diff options
Diffstat (limited to 'src/theory/arith/simplex-converge.cpp')
| -rw-r--r-- | src/theory/arith/simplex-converge.cpp | 1674 |
1 files changed, 0 insertions, 1674 deletions
diff --git a/src/theory/arith/simplex-converge.cpp b/src/theory/arith/simplex-converge.cpp deleted file mode 100644 index decce3882..000000000 --- a/src/theory/arith/simplex-converge.cpp +++ /dev/null @@ -1,1674 +0,0 @@ -/********************* */ -/*! \file simplex.cpp - ** \verbatim - ** Original author: taking - ** Major contributors: none - ** Minor contributors (to current version): none - ** This file is part of the CVC4 prototype. - ** Copyright (c) 2009, 2010, 2011 The Analysis of Computer Systems Group (ACSys) - ** Courant Institute of Mathematical Sciences - ** New York University - ** See the file COPYING in the top-level source directory for licensing - ** information.\endverbatim - ** - ** \brief [[ Add one-line brief description here ]] - ** - ** [[ Add lengthier description here ]] - ** \todo document this file - **/ - - -#include "theory/arith/simplex.h" -#include "theory/arith/options.h" - -using namespace std; - -using namespace CVC4; -using namespace CVC4::kind; - -using namespace CVC4::theory; -using namespace CVC4::theory::arith; - -static const bool CHECK_AFTER_PIVOT = true; - -SimplexDecisionProcedure::SimplexDecisionProcedure(LinearEqualityModule& linEq, NodeCallBack& conflictChannel, ArithVarMalloc& malloc, ConstraintDatabase& cd) : - d_conflictVariable(ARITHVAR_SENTINEL), - d_linEq(linEq), - d_partialModel(d_linEq.getPartialModel()), - d_tableau(d_linEq.getTableau()), - d_queue(d_partialModel, d_tableau), - d_numVariables(0), - d_conflictChannel(conflictChannel), - d_pivotsInRound(), - d_DELTA_ZERO(0,0), - d_arithVarMalloc(malloc), - d_constraintDatabase(cd), - d_optRow(ARITHVAR_SENTINEL), - d_negOptConstant(d_DELTA_ZERO) -{ - switch(ArithHeuristicPivotRule rule = options::arithHeuristicPivotRule()) { - case MINIMUM: - d_queue.setPivotRule(ArithPriorityQueue::MINIMUM); - break; - case BREAK_TIES: - d_queue.setPivotRule(ArithPriorityQueue::BREAK_TIES); - break; - case MAXIMUM: - d_queue.setPivotRule(ArithPriorityQueue::MAXIMUM); - break; - default: - Unhandled(rule); - } - - srand(62047); -} - -SimplexDecisionProcedure::Statistics::Statistics(): - d_statUpdateConflicts("theory::arith::UpdateConflicts", 0), - d_findConflictOnTheQueueTime("theory::arith::findConflictOnTheQueueTime"), - d_attemptBeforeDiffSearch("theory::arith::qi::BeforeDiffSearch::attempt",0), - d_successBeforeDiffSearch("theory::arith::qi::BeforeDiffSearch::success",0), - d_attemptAfterDiffSearch("theory::arith::qi::AfterDiffSearch::attempt",0), - d_successAfterDiffSearch("theory::arith::qi::AfterDiffSearch::success",0), - d_attemptDuringDiffSearch("theory::arith::qi::DuringDiffSearch::attempt",0), - d_successDuringDiffSearch("theory::arith::qi::DuringDiffSearch::success",0), - d_attemptDuringVarOrderSearch("theory::arith::qi::DuringVarOrderSearch::attempt",0), - d_successDuringVarOrderSearch("theory::arith::qi::DuringVarOrderSearch::success",0), - d_attemptAfterVarOrderSearch("theory::arith::qi::AfterVarOrderSearch::attempt",0), - d_successAfterVarOrderSearch("theory::arith::qi::AfterVarOrderSearch::success",0), - d_weakeningAttempts("theory::arith::weakening::attempts",0), - d_weakeningSuccesses("theory::arith::weakening::success",0), - d_weakenings("theory::arith::weakening::total",0), - d_weakenTime("theory::arith::weakening::time"), - d_simplexConflicts("theory::arith::simplexConflicts",0), - // primal - d_primalTimer("theory::arith::primal::overall::timer"), - d_internalTimer("theory::arith::primal::internal::timer"), - d_primalCalls("theory::arith::primal::calls",0), - d_primalSatCalls("theory::arith::primal::calls::sat",0), - d_primalUnsatCalls("theory::arith::primal::calls::unsat",0), - d_primalPivots("theory::arith::primal::pivots",0), - d_primalImprovingPivots("theory::arith::primal::pivots::improving",0), - d_primalThresholdReachedPivot("theory::arith::primal::thresholds",0), - d_primalThresholdReachedPivot_dropped("theory::arith::primal::thresholds::dropped",0), - d_primalReachedMaxPivots("theory::arith::primal::maxpivots",0), - d_primalReachedMaxPivots_contractMadeProgress("theory::arith::primal::maxpivots::contract",0), - d_primalReachedMaxPivots_checkForConflictWorked("theory::arith::primal::maxpivots::checkworked",0), - d_primalGlobalMinimum("theory::arith::primal::minimum",0), - d_primalGlobalMinimum_rowConflictWorked("theory::arith::primal::minimum::checkworked",0), - d_primalGlobalMinimum_firstHalfWasSat("theory::arith::primal::minimum::firsthalf::sat",0), - d_primalGlobalMinimum_firstHalfWasUnsat("theory::arith::primal::minimum::firsthalf::unsat",0), - d_primalGlobalMinimum_contractMadeProgress("theory::arith::primal::minimum::progress",0), - d_unboundedFound("theory::arith::primal::unbounded",0), - d_unboundedFound_drive("theory::arith::primal::unbounded::drive",0), - d_unboundedFound_dropped("theory::arith::primal::unbounded::dropped",0) -{ - StatisticsRegistry::registerStat(&d_statUpdateConflicts); - - StatisticsRegistry::registerStat(&d_findConflictOnTheQueueTime); - - StatisticsRegistry::registerStat(&d_attemptBeforeDiffSearch); - StatisticsRegistry::registerStat(&d_successBeforeDiffSearch); - StatisticsRegistry::registerStat(&d_attemptAfterDiffSearch); - StatisticsRegistry::registerStat(&d_successAfterDiffSearch); - StatisticsRegistry::registerStat(&d_attemptDuringDiffSearch); - StatisticsRegistry::registerStat(&d_successDuringDiffSearch); - StatisticsRegistry::registerStat(&d_attemptDuringVarOrderSearch); - StatisticsRegistry::registerStat(&d_successDuringVarOrderSearch); - StatisticsRegistry::registerStat(&d_attemptAfterVarOrderSearch); - StatisticsRegistry::registerStat(&d_successAfterVarOrderSearch); - - StatisticsRegistry::registerStat(&d_weakeningAttempts); - StatisticsRegistry::registerStat(&d_weakeningSuccesses); - StatisticsRegistry::registerStat(&d_weakenings); - StatisticsRegistry::registerStat(&d_weakenTime); - - StatisticsRegistry::registerStat(&d_simplexConflicts); - - //primal - StatisticsRegistry::registerStat(&d_primalTimer); - StatisticsRegistry::registerStat(&d_internalTimer); - - StatisticsRegistry::registerStat(&d_primalCalls); - StatisticsRegistry::registerStat(&d_primalSatCalls); - StatisticsRegistry::registerStat(&d_primalUnsatCalls); - - StatisticsRegistry::registerStat(&d_primalPivots); - StatisticsRegistry::registerStat(&d_primalImprovingPivots); - - StatisticsRegistry::registerStat(&d_primalThresholdReachedPivot); - StatisticsRegistry::registerStat(&d_primalThresholdReachedPivot_dropped); - - StatisticsRegistry::registerStat(&d_primalReachedMaxPivots); - StatisticsRegistry::registerStat(&d_primalReachedMaxPivots_contractMadeProgress); - StatisticsRegistry::registerStat(&d_primalReachedMaxPivots_checkForConflictWorked); - - - StatisticsRegistry::registerStat(&d_primalGlobalMinimum); - StatisticsRegistry::registerStat(&d_primalGlobalMinimum_rowConflictWorked); - StatisticsRegistry::registerStat(&d_primalGlobalMinimum_firstHalfWasSat); - StatisticsRegistry::registerStat(&d_primalGlobalMinimum_firstHalfWasUnsat); - StatisticsRegistry::registerStat(&d_primalGlobalMinimum_contractMadeProgress); - - - StatisticsRegistry::registerStat(&d_unboundedFound); - StatisticsRegistry::registerStat(&d_unboundedFound_drive); - StatisticsRegistry::registerStat(&d_unboundedFound_dropped); - -} - -SimplexDecisionProcedure::Statistics::~Statistics(){ - StatisticsRegistry::unregisterStat(&d_statUpdateConflicts); - - StatisticsRegistry::unregisterStat(&d_findConflictOnTheQueueTime); - - StatisticsRegistry::unregisterStat(&d_attemptBeforeDiffSearch); - StatisticsRegistry::unregisterStat(&d_successBeforeDiffSearch); - StatisticsRegistry::unregisterStat(&d_attemptAfterDiffSearch); - StatisticsRegistry::unregisterStat(&d_successAfterDiffSearch); - StatisticsRegistry::unregisterStat(&d_attemptDuringDiffSearch); - StatisticsRegistry::unregisterStat(&d_successDuringDiffSearch); - StatisticsRegistry::unregisterStat(&d_attemptDuringVarOrderSearch); - StatisticsRegistry::unregisterStat(&d_successDuringVarOrderSearch); - StatisticsRegistry::unregisterStat(&d_attemptAfterVarOrderSearch); - StatisticsRegistry::unregisterStat(&d_successAfterVarOrderSearch); - - StatisticsRegistry::unregisterStat(&d_weakeningAttempts); - StatisticsRegistry::unregisterStat(&d_weakeningSuccesses); - StatisticsRegistry::unregisterStat(&d_weakenings); - StatisticsRegistry::unregisterStat(&d_weakenTime); - - StatisticsRegistry::unregisterStat(&d_simplexConflicts); - - //primal - StatisticsRegistry::unregisterStat(&d_primalTimer); - StatisticsRegistry::unregisterStat(&d_internalTimer); - - StatisticsRegistry::unregisterStat(&d_primalCalls); - StatisticsRegistry::unregisterStat(&d_primalSatCalls); - StatisticsRegistry::unregisterStat(&d_primalUnsatCalls); - - StatisticsRegistry::unregisterStat(&d_primalPivots); - StatisticsRegistry::unregisterStat(&d_primalImprovingPivots); - - StatisticsRegistry::unregisterStat(&d_primalThresholdReachedPivot); - StatisticsRegistry::unregisterStat(&d_primalThresholdReachedPivot_dropped); - - StatisticsRegistry::unregisterStat(&d_primalReachedMaxPivots); - StatisticsRegistry::unregisterStat(&d_primalReachedMaxPivots_contractMadeProgress); - StatisticsRegistry::unregisterStat(&d_primalReachedMaxPivots_checkForConflictWorked); - - - StatisticsRegistry::unregisterStat(&d_primalGlobalMinimum); - StatisticsRegistry::unregisterStat(&d_primalGlobalMinimum_rowConflictWorked); - StatisticsRegistry::unregisterStat(&d_primalGlobalMinimum_firstHalfWasSat); - StatisticsRegistry::unregisterStat(&d_primalGlobalMinimum_firstHalfWasUnsat); - StatisticsRegistry::unregisterStat(&d_primalGlobalMinimum_contractMadeProgress); - - StatisticsRegistry::unregisterStat(&d_unboundedFound); - StatisticsRegistry::unregisterStat(&d_unboundedFound_drive); - StatisticsRegistry::unregisterStat(&d_unboundedFound_dropped); -} - - - - -ArithVar SimplexDecisionProcedure::minVarOrder(const SimplexDecisionProcedure& simp, ArithVar x, ArithVar y){ - Assert(x != ARITHVAR_SENTINEL); - Assert(y != ARITHVAR_SENTINEL); - // Assert(!simp.d_tableau.isBasic(x)); - // Assert(!simp.d_tableau.isBasic(y)); - if(x <= y){ - return x; - } else { - return y; - } -} - -ArithVar SimplexDecisionProcedure::minColLength(const SimplexDecisionProcedure& simp, ArithVar x, ArithVar y){ - Assert(x != ARITHVAR_SENTINEL); - Assert(y != ARITHVAR_SENTINEL); - Assert(!simp.d_tableau.isBasic(x)); - Assert(!simp.d_tableau.isBasic(y)); - uint32_t xLen = simp.d_tableau.getColLength(x); - uint32_t yLen = simp.d_tableau.getColLength(y); - if( xLen > yLen){ - return y; - } else if( xLen== yLen ){ - return minVarOrder(simp,x,y); - }else{ - return x; - } -} - -ArithVar SimplexDecisionProcedure::minRowLength(const SimplexDecisionProcedure& simp, ArithVar x, ArithVar y){ - Assert(x != ARITHVAR_SENTINEL); - Assert(y != ARITHVAR_SENTINEL); - Assert(simp.d_tableau.isBasic(x)); - Assert(simp.d_tableau.isBasic(y)); - uint32_t xLen = simp.d_tableau.getRowLength(simp.d_tableau.basicToRowIndex(x)); - uint32_t yLen = simp.d_tableau.getRowLength(simp.d_tableau.basicToRowIndex(y)); - if( xLen > yLen){ - return y; - } else if( xLen == yLen ){ - return minVarOrder(simp,x,y); - }else{ - return x; - } -} -ArithVar SimplexDecisionProcedure::minBoundAndRowCount(const SimplexDecisionProcedure& simp, ArithVar x, ArithVar y){ - Assert(x != ARITHVAR_SENTINEL); - Assert(y != ARITHVAR_SENTINEL); - Assert(!simp.d_tableau.isBasic(x)); - Assert(!simp.d_tableau.isBasic(y)); - if(simp.d_partialModel.hasEitherBound(x) && !simp.d_partialModel.hasEitherBound(y)){ - return y; - }else if(!simp.d_partialModel.hasEitherBound(x) && simp.d_partialModel.hasEitherBound(y)){ - return x; - }else { - return minColLength(simp, x, y); - } -} - -template <bool above> -ArithVar SimplexDecisionProcedure::selectSlack(ArithVar x_i, SimplexDecisionProcedure::PreferenceFunction pref){ - ArithVar slack = ARITHVAR_SENTINEL; - - for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd(); ++iter){ - const Tableau::Entry& entry = *iter; - ArithVar nonbasic = entry.getColVar(); - if(nonbasic == x_i) continue; - - const Rational& a_ij = entry.getCoefficient(); - int sgn = a_ij.sgn(); - if(isAcceptableSlack<above>(sgn, nonbasic)){ - //If one of the above conditions is met, we have found an acceptable - //nonbasic variable to pivot x_i with. We can now choose which one we - //prefer the most. - slack = (slack == ARITHVAR_SENTINEL) ? nonbasic : pref(*this, slack, nonbasic); - } - } - - return slack; -} - -Node betterConflict(TNode x, TNode y){ - if(x.isNull()) return y; - else if(y.isNull()) return x; - else if(x.getNumChildren() <= y.getNumChildren()) return x; - else return y; -} - -bool SimplexDecisionProcedure::findConflictOnTheQueue(SearchPeriod type) { - TimerStat::CodeTimer codeTimer(d_statistics.d_findConflictOnTheQueueTime); - Assert(d_successes.empty()); - - switch(type){ - case BeforeDiffSearch: ++(d_statistics.d_attemptBeforeDiffSearch); break; - case DuringDiffSearch: ++(d_statistics.d_attemptDuringDiffSearch); break; - case AfterDiffSearch: ++(d_statistics.d_attemptAfterDiffSearch); break; - case DuringVarOrderSearch: ++(d_statistics.d_attemptDuringVarOrderSearch); break; - case AfterVarOrderSearch: ++(d_statistics.d_attemptAfterVarOrderSearch); break; - } - - ArithPriorityQueue::const_iterator i = d_queue.begin(); - ArithPriorityQueue::const_iterator end = d_queue.end(); - for(; i != end; ++i){ - ArithVar x_i = *i; - - if(x_i != d_conflictVariable && d_tableau.isBasic(x_i) && !d_successes.isMember(x_i)){ - Node possibleConflict = checkBasicForConflict(x_i); - if(!possibleConflict.isNull()){ - d_successes.add(x_i); - reportConflict(possibleConflict); - } - } - } - if(!d_successes.empty()){ - switch(type){ - case BeforeDiffSearch: ++(d_statistics.d_successBeforeDiffSearch); break; - case DuringDiffSearch: ++(d_statistics.d_successDuringDiffSearch); break; - case AfterDiffSearch: ++(d_statistics.d_successAfterDiffSearch); break; - case DuringVarOrderSearch: ++(d_statistics.d_successDuringVarOrderSearch); break; - case AfterVarOrderSearch: ++(d_statistics.d_successAfterVarOrderSearch); break; - } - d_successes.purge(); - return true; - }else{ - return false; - } -} - -Result::Sat SimplexDecisionProcedure::dualFindModel(bool exactResult){ - Assert(d_conflictVariable == ARITHVAR_SENTINEL); - Assert(d_queue.inCollectionMode()); - - if(d_queue.empty()){ - return Result::SAT; - } - static CVC4_THREADLOCAL(unsigned int) instance = 0; - instance = instance + 1; - Debug("arith::findModel") << "begin findModel()" << instance << endl; - - d_queue.transitionToDifferenceMode(); - - Result::Sat result = Result::SAT_UNKNOWN; - - if(d_queue.empty()){ - result = Result::SAT; - }else if(d_queue.size() > 1){ - if(findConflictOnTheQueue(BeforeDiffSearch)){ - result = Result::UNSAT; - } - } - static const bool verbose = false; - exactResult |= options::arithStandardCheckVarOrderPivots() < 0; - const uint32_t inexactResultsVarOrderPivots = exactResult ? 0 : options::arithStandardCheckVarOrderPivots(); - - uint32_t checkPeriod = options::arithSimplexCheckPeriod(); - if(result == Result::SAT_UNKNOWN){ - uint32_t numDifferencePivots = options::arithHeuristicPivots() < 0 ? - d_numVariables + 1 : options::arithHeuristicPivots(); - // The signed to unsigned conversion is safe. - uint32_t pivotsRemaining = numDifferencePivots; - while(!d_queue.empty() && - result != Result::UNSAT && - pivotsRemaining > 0){ - uint32_t pivotsToDo = min(checkPeriod, pivotsRemaining); - pivotsRemaining -= pivotsToDo; - if(searchForFeasibleSolution(pivotsToDo)){ - result = Result::UNSAT; - }//Once every CHECK_PERIOD examine the entire queue for conflicts - if(result != Result::UNSAT){ - if(findConflictOnTheQueue(DuringDiffSearch)) { result = Result::UNSAT; } - }else{ - findConflictOnTheQueue(AfterDiffSearch); // already unsat - } - } - - if(verbose && numDifferencePivots > 0){ - if(result == Result::UNSAT){ - Message() << "diff order found unsat" << endl; - }else if(d_queue.empty()){ - Message() << "diff order found model" << endl; - }else{ - Message() << "diff order missed" << endl; - } - } - } - - if(!d_queue.empty() && result != Result::UNSAT){ - if(exactResult){ - d_queue.transitionToVariableOrderMode(); - - while(!d_queue.empty() && result != Result::UNSAT){ - if(searchForFeasibleSolution(checkPeriod)){ - result = Result::UNSAT; - } - - //Once every CHECK_PERIOD examine the entire queue for conflicts - if(result != Result::UNSAT){ - if(findConflictOnTheQueue(DuringVarOrderSearch)){ - result = Result::UNSAT; - } - } else{ - findConflictOnTheQueue(AfterVarOrderSearch); - } - } - if(verbose){ - if(result == Result::UNSAT){ - Message() << "bland found unsat" << endl; - }else if(d_queue.empty()){ - Message() << "bland found model" << endl; - }else{ - Message() << "bland order missed" << endl; - } - } - }else{ - d_queue.transitionToVariableOrderMode(); - - if(searchForFeasibleSolution(inexactResultsVarOrderPivots)){ - result = Result::UNSAT; - findConflictOnTheQueue(AfterVarOrderSearch); // already unsat - }else{ - if(findConflictOnTheQueue(AfterVarOrderSearch)) { result = Result::UNSAT; } - } - - if(verbose){ - if(result == Result::UNSAT){ - Message() << "restricted var order found unsat" << endl; - }else if(d_queue.empty()){ - Message() << "restricted var order found model" << endl; - }else{ - Message() << "restricted var order missed" << endl; - } - } - } - } - - if(result == Result::SAT_UNKNOWN && d_queue.empty()){ - result = Result::SAT; - } - - - - d_pivotsInRound.purge(); - // ensure that the conflict variable is still in the queue. - if(d_conflictVariable != ARITHVAR_SENTINEL){ - d_queue.enqueueIfInconsistent(d_conflictVariable); - } - d_conflictVariable = ARITHVAR_SENTINEL; - - d_queue.transitionToCollectionMode(); - Assert(d_queue.inCollectionMode()); - Debug("arith::findModel") << "end findModel() " << instance << " " << result << endl; - return result; - - - // Assert(foundConflict || d_queue.empty()); - - // // Curiously the invariant that we always do a full check - // // means that the assignment we can always empty these queues. - // d_queue.clear(); - // d_pivotsInRound.purge(); - // d_conflictVariable = ARITHVAR_SENTINEL; - - // Assert(!d_queue.inCollectionMode()); - // d_queue.transitionToCollectionMode(); - - - // Assert(d_queue.inCollectionMode()); - - // Debug("arith::findModel") << "end findModel() " << instance << endl; - - // return foundConflict; -} - - - -Node SimplexDecisionProcedure::checkBasicForConflict(ArithVar basic){ - - Assert(d_tableau.isBasic(basic)); - const DeltaRational& beta = d_partialModel.getAssignment(basic); - - if(d_partialModel.strictlyLessThanLowerBound(basic, beta)){ - ArithVar x_j = selectSlackUpperBound(basic); - if(x_j == ARITHVAR_SENTINEL ){ - return generateConflictBelowLowerBound(basic); - } - }else if(d_partialModel.strictlyGreaterThanUpperBound(basic, beta)){ - ArithVar x_j = selectSlackLowerBound(basic); - if(x_j == ARITHVAR_SENTINEL ){ - return generateConflictAboveUpperBound(basic); - } - } - return Node::null(); -} - -//corresponds to Check() in dM06 -//template <SimplexDecisionProcedure::PreferenceFunction pf> -bool SimplexDecisionProcedure::searchForFeasibleSolution(uint32_t remainingIterations){ - Debug("arith") << "searchForFeasibleSolution" << endl; - Assert(remainingIterations > 0); - - while(remainingIterations > 0){ - if(Debug.isOn("paranoid:check_tableau")){ d_linEq.debugCheckTableau(); } - - ArithVar x_i = d_queue.dequeueInconsistentBasicVariable(); - Debug("arith::update::select") << "selectSmallestInconsistentVar()=" << x_i << endl; - if(x_i == ARITHVAR_SENTINEL){ - Debug("arith_update") << "No inconsistent variables" << endl; - return false; //sat - } - - --remainingIterations; - - bool useVarOrderPivot = d_pivotsInRound.count(x_i) >= options::arithPivotThreshold(); - if(!useVarOrderPivot){ - d_pivotsInRound.add(x_i); - } - - - Debug("playground") << "pivots in rounds: " << d_pivotsInRound.count(x_i) - << " use " << useVarOrderPivot - << " threshold " << options::arithPivotThreshold() - << endl; - - PreferenceFunction pf = useVarOrderPivot ? minVarOrder : minBoundAndRowCount; - - DeltaRational beta_i = d_partialModel.getAssignment(x_i); - ArithVar x_j = ARITHVAR_SENTINEL; - - if(d_partialModel.strictlyLessThanLowerBound(x_i, beta_i)){ - x_j = selectSlackUpperBound(x_i, pf); - if(x_j == ARITHVAR_SENTINEL ){ - ++(d_statistics.d_statUpdateConflicts); - Node conflict = generateConflictBelowLowerBound(x_i); //unsat - d_conflictVariable = x_i; - reportConflict(conflict); - return true; - } - DeltaRational l_i = d_partialModel.getLowerBound(x_i); - d_linEq.pivotAndUpdate(x_i, x_j, l_i); - - }else if(d_partialModel.strictlyGreaterThanUpperBound(x_i, beta_i)){ - x_j = selectSlackLowerBound(x_i, pf); - if(x_j == ARITHVAR_SENTINEL ){ - ++(d_statistics.d_statUpdateConflicts); - Node conflict = generateConflictAboveUpperBound(x_i); //unsat - - d_conflictVariable = x_i; - reportConflict(conflict); - return true; - } - DeltaRational u_i = d_partialModel.getUpperBound(x_i); - d_linEq.pivotAndUpdate(x_i, x_j, u_i); - } - Assert(x_j != ARITHVAR_SENTINEL); - - //Check to see if we already have a conflict with x_j to prevent wasteful work - if(CHECK_AFTER_PIVOT){ - Node possibleConflict = checkBasicForConflict(x_j); - if(!possibleConflict.isNull()){ - d_conflictVariable = x_j; - reportConflict(possibleConflict); - return true; // unsat - } - } - } - Assert(remainingIterations == 0); - - return false; -} - - - -Constraint SimplexDecisionProcedure::weakestExplanation(bool aboveUpper, DeltaRational& surplus, ArithVar v, const Rational& coeff, bool& anyWeakening, ArithVar basic){ - - int sgn = coeff.sgn(); - bool ub = aboveUpper?(sgn < 0) : (sgn > 0); - - Constraint c = ub ? - d_partialModel.getUpperBoundConstraint(v) : - d_partialModel.getLowerBoundConstraint(v); - -// #warning "revisit" -// Node exp = ub ? -// d_partialModel.explainUpperBound(v) : -// d_partialModel.explainLowerBound(v); - - bool weakened; - do{ - const DeltaRational& bound = c->getValue(); - - weakened = false; - - Constraint weaker = ub? - c->getStrictlyWeakerUpperBound(true, true): - c->getStrictlyWeakerLowerBound(true, true); - - // Node weaker = ub? - // d_propManager.strictlyWeakerAssertedUpperBound(v, bound): - // d_propManager.strictlyWeakerAssertedLowerBound(v, bound); - - if(weaker != NullConstraint){ - //if(!weaker.isNull()){ - const DeltaRational& weakerBound = weaker->getValue(); - //DeltaRational weakerBound = asDeltaRational(weaker); - - DeltaRational diff = aboveUpper ? bound - weakerBound : weakerBound - bound; - //if var == basic, - // if aboveUpper, weakerBound > bound, multiply by -1 - // if !aboveUpper, weakerBound < bound, multiply by -1 - diff = diff * coeff; - if(surplus > diff){ - ++d_statistics.d_weakenings; - weakened = true; - anyWeakening = true; - surplus = surplus - diff; - - Debug("weak") << "found:" << endl; - if(v == basic){ - Debug("weak") << " basic: "; - } - Debug("weak") << " " << surplus << " "<< diff << endl - << " " << bound << c << endl - << " " << weakerBound << weaker << endl; - - Assert(diff > d_DELTA_ZERO); - c = weaker; - } - } - }while(weakened); - - return c; -} - -Node SimplexDecisionProcedure::weakenConflict(bool aboveUpper, ArithVar basicVar){ - TimerStat::CodeTimer codeTimer(d_statistics.d_weakenTime); - - const DeltaRational& assignment = d_partialModel.getAssignment(basicVar); - DeltaRational surplus; - if(aboveUpper){ - Assert(d_partialModel.hasUpperBound(basicVar)); - Assert(assignment > d_partialModel.getUpperBound(basicVar)); - surplus = assignment - d_partialModel.getUpperBound(basicVar); - }else{ - Assert(d_partialModel.hasLowerBound(basicVar)); - Assert(assignment < d_partialModel.getLowerBound(basicVar)); - surplus = d_partialModel.getLowerBound(basicVar) - assignment; - } - - NodeBuilder<> conflict(kind::AND); - bool anyWeakenings = false; - for(Tableau::RowIterator i = d_tableau.basicRowIterator(basicVar); !i.atEnd(); ++i){ - const Tableau::Entry& entry = *i; - ArithVar v = entry.getColVar(); - const Rational& coeff = entry.getCoefficient(); - bool weakening = false; - Constraint c = weakestExplanation(aboveUpper, surplus, v, coeff, weakening, basicVar); - Debug("weak") << "weak : " << weakening << " " - << c->assertedToTheTheory() << " " - << d_partialModel.getAssignment(v) << " " - << c << endl - << c->explainForConflict() << endl; - anyWeakenings = anyWeakenings || weakening; - - Debug("weak") << "weak : " << c->explainForConflict() << endl; - c->explainForConflict(conflict); - } - ++d_statistics.d_weakeningAttempts; - if(anyWeakenings){ - ++d_statistics.d_weakeningSuccesses; - } - return conflict; -} - - -Node SimplexDecisionProcedure::generateConflictAboveUpperBound(ArithVar conflictVar){ - return weakenConflict(true, conflictVar); -} - -Node SimplexDecisionProcedure::generateConflictBelowLowerBound(ArithVar conflictVar){ - return weakenConflict(false, conflictVar); -} - - -// responses -// unbounded below(arithvar) -// reached threshold -// reached maxpivots -// reached GlobalMinimumd -// -SimplexDecisionProcedure::PrimalResponse SimplexDecisionProcedure::primal(uint32_t maxIterations) -{ - Assert(d_optRow != ARITHVAR_SENTINEL); - - for(uint32_t iteration = 0; iteration < maxIterations; iteration++){ - if(belowThreshold()){ - return ReachedThresholdValue; - } - - PrimalResponse res = primalCheck(); - switch(res){ - case GlobalMinimum: - case FoundUnboundedVariable: - return res; - case NoProgressOnLeaving: - ++d_statistics.d_primalPivots; - ++d_pivotsSinceOptProgress; - ++d_pivotsSinceLastCheck; - ++d_pivotsSinceErrorProgress; - - d_linEq.pivotAndUpdate(d_primalCarry.d_entering, d_primalCarry.d_leaving, d_partialModel.getAssignment(d_primalCarry.d_entering)); - - if(Debug.isOn("primal::tableau")){ - d_linEq.debugCheckTableau(); - } - if(Debug.isOn("primal::consistent")){ Assert(d_linEq.debugEntireLinEqIsConsistent("MakeProgressOnLeaving")); } - - break; - case MakeProgressOnLeaving: - { - ++d_statistics.d_primalPivots; - ++d_statistics.d_primalImprovingPivots; - - d_pivotsSinceOptProgress = 0; - ++d_pivotsSinceErrorProgress; - ++d_pivotsSinceLastCheck; - - d_linEq.pivotAndUpdate(d_primalCarry.d_entering, d_primalCarry.d_leaving, d_primalCarry.d_nextEnteringValue); - - static int helpful = 0; - static int hurtful = 0; - static int penguin = 0; - if(d_currentErrorVariables.isKey(d_primalCarry.d_entering)){ - cout << "entering is error " << d_primalCarry.d_entering; - if(d_currentErrorVariables[d_primalCarry.d_entering].errorIsLeqZero(d_partialModel)){ - cout << " now below threshold (" << (++helpful) << ") " << d_pivotsSinceErrorProgress << endl; - }else{ - cout << "ouch (" << (++hurtful) << ")"<< d_pivotsSinceErrorProgress << endl; - } - }else if(d_currentErrorVariables.isKey(d_primalCarry.d_leaving)){ - cout << "leaving is error " << d_primalCarry.d_leaving; - if(d_currentErrorVariables[d_primalCarry.d_leaving].errorIsLeqZero(d_partialModel)){ - cout << " now below threshold(" << (++helpful) << ")" << d_pivotsSinceErrorProgress << endl; - }else{ - cout << "ouch (" << (++hurtful) << ")" << d_pivotsSinceErrorProgress<< endl; - } - }else{ - cout << " penguin (" << (++penguin) << ")" << d_pivotsSinceErrorProgress<< endl; - } - - if(Debug.isOn("primal::tableau")){ - d_linEq.debugCheckTableau(); - } - if(Debug.isOn("primal::consistent")){ Assert(d_linEq.debugEntireLinEqIsConsistent("MakeProgressOnLeaving")); } - } - break; - default: - Unreachable(); - } - } - return UsedMaxPivots; -} - - -/** - * Error set - * ErrorVariable |-> {ErrorVariable, InputVariable, InputConstraint} - */ - -/** - * Returns SAT if it was able to satisfy all of the constraints in the error set - * Returns UNSAT if it was able to able to find an error - */ -Result::Sat SimplexDecisionProcedure::primalConverge(int depth){ - d_pivotsSinceLastCheck = 0; - - while(!d_currentErrorVariables.empty()){ - PrimalResponse res = primal(MAX_ITERATIONS - d_pivotsSinceLastCheck); - - switch(res){ - case FoundUnboundedVariable: - - // Drive the variable to at least 0 - // TODO This variable should be driven to a value that makes all of the error functions including it 0 - // It'll or another unbounded will be selected in the next round anyways so ignore for now. - ++d_statistics.d_unboundedFound; - if( !belowThreshold() ){ - driveOptToZero(d_primalCarry.d_unbounded); - d_linEq.debugCheckTableau(); - if(Debug.isOn("primal::consistent")){ Assert(d_linEq.debugEntireLinEqIsConsistent("primalConverge")); } - - ++d_statistics.d_unboundedFound_drive; - } - Assert(belowThreshold()); - { - uint32_t dropped = contractErrorVariables(true); - Debug("primal::converge") << "primalConverge -> FoundUnboundedVariable -> dropped " << dropped << " to " << d_currentErrorVariables.size() << endl; - d_statistics.d_unboundedFound_dropped += dropped; - } - break; - case ReachedThresholdValue: - ++d_statistics.d_primalThresholdReachedPivot; - - Assert(belowThreshold()); - { - uint32_t dropped = contractErrorVariables(true); - Debug("primal::converge") << "primalConverge -> ReachedThresholdValue -> dropped " << dropped << " to " << d_currentErrorVariables.size() << endl; - d_statistics.d_primalThresholdReachedPivot_dropped += dropped; - } - break; - case UsedMaxPivots: - { - d_pivotsSinceLastCheck = 0; - ++d_statistics.d_primalReachedMaxPivots; - - // periodically attempt to do the following : - // contract the error variable - // check for a conflict on an error variable - uint32_t dropped = contractErrorVariables(false); - - if( checkForRowConflicts() ){ // try to periodically look for a row - Debug("primal::converge") << "primalConverge -> UsedMaxPivots -> unsat " << endl; - - ++d_statistics.d_primalReachedMaxPivots_checkForConflictWorked; - return Result::UNSAT; // row conflicts are minimal so stop. - } - - if(dropped > 0){ - Debug("primal::converge") << "primalConverge -> UsedMaxPivots -> dropped " << dropped << " to " << d_currentErrorVariables.size() << endl; - ++d_statistics.d_primalReachedMaxPivots_contractMadeProgress; - }else{ - Debug("primal::converge") << "primalConverge -> UsedMaxPivots -> nothing " << endl; - } - } - break; - case GlobalMinimum: - ++d_statistics.d_primalGlobalMinimum; - - // If the minimum is positive, this is unsat. - // However, the optimization row is not necessarily a minimal conflict - if(!belowThreshold()){ - - if(d_currentErrorVariables.size() == 1){ - // The optimization function is exactly the same as the last remaining variable - // The conflict for the row is the same as the conflict for the optimization function. - bool foundConflict = checkForRowConflicts(); - Assert(foundConflict); - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> one variable" << endl; - - return Result::UNSAT; - }else{ - // There are at least 2 error variables. - // Look for a minimal conflict - - - if(checkForRowConflicts() ){ - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> postitive -> rowconflict " << endl; - - ++d_statistics.d_primalGlobalMinimum_rowConflictWorked; - return Result::UNSAT; - } - - uint32_t dropped = contractErrorVariables(false); - - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> postitive -> dropped " << dropped << " to " << d_currentErrorVariables.size() << endl; - if(dropped > 0){ - ++d_statistics.d_primalGlobalMinimum_contractMadeProgress; - } - - ErrorMap half; - d_currentErrorVariables.splitInto(half); - - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> recursion " << depth << endl; - - - reconstructOptimizationFunction(); - Result::Sat resultOnRemaining = primalConverge(depth + 1); - - if(resultOnRemaining == Result::UNSAT){ - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> recursion " << depth << " was unsat " << endl; - ++d_statistics.d_primalGlobalMinimum_firstHalfWasUnsat; - restoreErrorVariables(half); - return Result::UNSAT; - }else{ - ++d_statistics.d_primalGlobalMinimum_firstHalfWasSat; - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> recursion " << depth << " was sat " << endl; - - Assert(resultOnRemaining == Result::SAT); - Assert(d_currentErrorVariables.empty()); - d_currentErrorVariables.addAll(half); - reconstructOptimizationFunction(); - return primalConverge(depth + 1); - } - } - - }else{ - - // if the optimum is <= 0 - // drop all of the satisfied variables and continue; - uint32_t dropped = contractErrorVariables(true); - Debug("primal::converge") << "primalConverge -> GlobalMinimum -> negative -> dropped "<< dropped << " to " << d_currentErrorVariables.size() << endl; - - ++d_statistics.d_primalGlobalMinimum_contractMadeProgress; - } - break; - default: - Unreachable(); - } - } - - return Result::SAT; -} - - -Result::Sat SimplexDecisionProcedure::primalFindModel(){ - Assert(d_primalCarry.isClear()); - - // Reduce the queue to only contain violations - reduceQueue(); - - if(d_queue.empty()){ - return Result::SAT; - } - TimerStat::CodeTimer codeTimer(d_statistics.d_primalTimer); - - ++d_statistics.d_primalCalls; - - Debug("primalFindModel") << "primalFindModel() begin" << endl; - - const int PAUSE_RATE = 100; - if(Debug.isOn("primal::pause") && d_statistics.d_primalCalls.getData() % PAUSE_RATE == 0){ - Debug("primal::pause") << "waiting for input: "; - std::string dummy; - std::getline(std::cin, dummy); - } - - // TODO restore the tableau by ejecting variables - // Tableau copy(d_tableau); - - Result::Sat answer; - { - TimerStat::CodeTimer codeTimer(d_statistics.d_internalTimer); - - // This is needed because of the fiddling with the partial model - //context::Context::ScopedPush speculativePush(satContext); - - constructErrorVariables(); - constructOptimizationFunction(); - if(Debug.isOn("primal::tableau")){ d_linEq.debugCheckTableau(); } - if(Debug.isOn("primal::consistent")){ d_linEq.debugEntireLinEqIsConsistent("primalFindModel 1"); } - answer = primalConverge(0); - } - removeOptimizationFunction(); - - - // exit - uint32_t nc = d_tableau.getNumColumns(); - //d_tableau = copy; - while(d_tableau.getNumColumns() < nc){ - d_tableau.increaseSize(); - } - restoreErrorVariables(d_currentErrorVariables); - - reduceQueue(); - - if(Debug.isOn("primal::tableau")){ d_linEq.debugCheckTableau(); } - - if(Debug.isOn("primal::consistent")){ d_linEq.debugEntireLinEqIsConsistent("primalFindModel2"); } - Debug("primalFindModel") << "primalFindModel() end " << answer << endl; - - // The set of variables in conflict with their bounds will still be a subset of the - // variables that are in conflict with their bounds in the beginning. - // The basic variables are the same because of the copy. - // Thus it is safe to not try to not recompute the queue of violating variables - - if(answer == Result::UNSAT){ - // This needs to be done in a different context level than the push - ++d_statistics.d_primalUnsatCalls; - }else{ - ++d_statistics.d_primalSatCalls; - } - - d_primalCarry.clear(); - - return answer; -} - -/** Clears the ErrorMap and relase the resources associated with it. - * There are a couple of error maps around - */ -void SimplexDecisionProcedure::restoreErrorVariables(SimplexDecisionProcedure::ErrorMap& es){ - while(!es.empty()){ - ArithVar e = es.back(); - - reassertErrorConstraint(e, es, false, false); - es.pop_back(); - } -} - -void SimplexDecisionProcedure::constructErrorVariables(){ - Assert(d_currentErrorVariables.empty()); - Assert(!d_queue.empty()); - - for(ArithPriorityQueue::const_iterator iter = d_queue.begin(), end = d_queue.end(); iter != end; ++iter){ - ArithVar input = *iter; - - Assert(d_tableau.isBasic(input)); - Assert(!d_partialModel.assignmentIsConsistent(input)); - - Assert(!d_currentErrorVariables.isKey(input)); - - bool ub = d_partialModel.strictlyGreaterThanUpperBound(input, d_partialModel.getAssignment(input)); - - Constraint original = ub ? d_partialModel.getUpperBoundConstraint(input) - : d_partialModel.getLowerBoundConstraint(input); - - d_currentErrorVariables.set(input, ErrorInfo(input, ub, original)); - - if(ub){ - d_partialModel.forceRelaxUpperBound(input); - }else{ - d_partialModel.forceRelaxLowerBound(input); - } - Debug("primal") << "adding error variable (" << input << ", " << ", " << original <<") "; - Debug("primal") << "ub "<< ub << " " << d_partialModel.getAssignment(input) << " " << original->getValue() <<")" << endl; - d_currentErrorVariables.set(input, ErrorInfo(input, ub, original)); - - // Constraint boundIsValue = d_constraintDatabase.getConstraint(bound, Equality, original->getValue()); - // boundIsValue->setPsuedoConstraint(); - - // d_partialModel.setAssignment(bound, boundIsValue->getValue()); - // d_partialModel.setUpperBoundConstraint(boundIsValue); - // d_partialModel.setLowerBoundConstraint(boundIsValue); - - // // if ub - // // then error = x - boundIsValue - // // else error = boundIsValue - x - - // ArithVar error = requestVariable(); - - // DeltaRational diff = ub ? - // d_partialModel.getAssignment(input) - boundIsValue->getValue() : - // boundIsValue->getValue() - d_partialModel.getAssignment(input); - - // d_partialModel.setAssignment(error, diff); - - // vector<Rational> coeffs; - // vector<ArithVar> variables; - // variables.push_back(input); - // coeffs.push_back(ub ? Rational(1) : Rational(-1)); - // variables.push_back(bound); - // coeffs.push_back(ub ? Rational(-1) : Rational(1)); - - // d_tableau.addRow(error, coeffs, variables); - - - } - - if(Debug.isOn("primal::tableau")){ d_linEq.debugCheckTableau(); } - if(Debug.isOn("primal::consistent")){ d_linEq.debugEntireLinEqIsConsistent("constructErrorVariables");} - Assert(!d_currentErrorVariables.empty()); -} - - - -/** Returns true if it has found a row conflict for any of the error variables. */ -bool SimplexDecisionProcedure::checkForRowConflicts(){ - vector<ArithVar> inConflict; - for(ErrorMap::const_iterator iter = d_currentErrorVariables.begin(), end = d_currentErrorVariables.end(); iter != end; ++iter){ - ArithVar error = *iter; - const ErrorInfo& info = d_currentErrorVariables[error]; - if(d_tableau.isBasic(error) && !info.errorIsLeqZero(d_partialModel)){ - - ArithVar x_j = info.isUpperbound() ? - selectSlackLowerBound(error) : - selectSlackUpperBound(error); - - if(x_j == ARITHVAR_SENTINEL ){ - inConflict.push_back(error); - } - } - } - - if(!inConflict.empty()){ - while(!inConflict.empty()){ - ArithVar error = inConflict.back(); - inConflict.pop_back(); - - reassertErrorConstraint(error, d_currentErrorVariables, false, true); - - Node conflict = d_currentErrorVariables[error].isUpperbound() ? - generateConflictAboveUpperBound(error) : - generateConflictBelowLowerBound(error); - Assert(!conflict.isNull()); - - d_currentErrorVariables.remove(error); - - reportConflict(conflict); - } - reconstructOptimizationFunction(); - return true; - }else{ - return false; - } -} - -void SimplexDecisionProcedure::reassertErrorConstraint(ArithVar e, SimplexDecisionProcedure::ErrorMap& es, bool definitelyTrue, bool definitelyFalse){ - Assert(es.isKey(e)); - const ErrorInfo& info = es.get(e); - Constraint original = info.getConstraint(); - - if(info.isUpperbound()){ - d_partialModel.setUpperBoundConstraint(original); - }else if(original->isLowerBound()){ - d_partialModel.setLowerBoundConstraint(original); - } - - Assert(!definitelyTrue || d_partialModel.assignmentIsConsistent(e)); - Assert(!definitelyFalse || !d_partialModel.assignmentIsConsistent(e)); -} - -uint32_t SimplexDecisionProcedure::contractErrorVariables(bool guaranteedSuccess){ - uint32_t entrySize = d_currentErrorVariables.size(); - Debug("primal::contract") << "contractErrorVariables() begin : " << d_currentErrorVariables.size() << endl; - - std::vector<ArithVar> toRemove; - for(ErrorMap::const_iterator iter = d_currentErrorVariables.begin(), end = d_currentErrorVariables.end(); iter != end; ++iter){ - ArithVar e = *iter; - if(d_currentErrorVariables[e].errorIsLeqZero(d_partialModel)){ - toRemove.push_back(e); - } - } - - Assert(!guaranteedSuccess || !toRemove.empty()); - - if(!toRemove.empty()){ - while(!toRemove.empty()){ - ArithVar e = toRemove.back(); - toRemove.pop_back(); - reassertErrorConstraint(e, d_currentErrorVariables, true, false); - d_currentErrorVariables.remove(e); - } - - reconstructOptimizationFunction(); - } - - Debug("primal::contract") << "contractErrorVariables() end : " << d_currentErrorVariables.size() << endl; - - uint32_t exitSize = d_currentErrorVariables.size(); - - Assert(exitSize <= entrySize); - Assert(!guaranteedSuccess|| exitSize < entrySize); - return entrySize - exitSize; -} - -void SimplexDecisionProcedure::removeOptimizationFunction(){ - Assert(d_optRow != ARITHVAR_SENTINEL); - Assert(d_tableau.isBasic(d_optRow)); - - d_tableau.removeBasicRow(d_optRow); - releaseVariable(d_optRow); - - d_optRow = ARITHVAR_SENTINEL; - d_negOptConstant = d_DELTA_ZERO; - - Assert(d_optRow == ARITHVAR_SENTINEL); -} - -void SimplexDecisionProcedure::constructOptimizationFunction(){ - Assert(d_optRow == ARITHVAR_SENTINEL); - Assert(d_negOptConstant == d_DELTA_ZERO); - - d_optRow = requestVariable(); - - - std::vector<Rational> coeffs; - std::vector<ArithVar> variables; - for(ErrorMap::const_iterator iter = d_currentErrorVariables.begin(), end = d_currentErrorVariables.end(); iter != end; ++iter){ - ArithVar e = *iter; - - if(d_currentErrorVariables[e].isUpperbound()){ - coeffs.push_back(Rational(1)); - variables.push_back(e); - d_negOptConstant = d_negOptConstant + (d_currentErrorVariables[e].getConstraint()->getValue()); - }else{ - coeffs.push_back(Rational(-1)); - variables.push_back(e); - d_negOptConstant = d_negOptConstant - (d_currentErrorVariables[e].getConstraint()->getValue()); - } - } - d_tableau.addRow(d_optRow, coeffs, variables); - - DeltaRational newAssignment = d_linEq.computeRowValue(d_optRow, false); - d_partialModel.setAssignment(d_optRow, newAssignment); - - if(Debug.isOn("primal::tableau")){ d_linEq.debugCheckTableau(); } - - - if(Debug.isOn("primal::consistent")){ - d_linEq.debugEntireLinEqIsConsistent("constructOptimizationFunction"); - } - - d_pivotsSinceOptProgress = 0; - d_pivotsSinceErrorProgress = 0; - - Assert(d_optRow != ARITHVAR_SENTINEL); -} - -void SimplexDecisionProcedure::reconstructOptimizationFunction(){ - removeOptimizationFunction(); - constructOptimizationFunction(); -} - - - -/* TODO: - * Very naive implementation. Recomputes everything every time. - * Currently looks for the variable that can decrease the optimization function the most. - * - */ -SimplexDecisionProcedure::PrimalResponse SimplexDecisionProcedure::primalCheck() -{ - Debug("primal") << "primalCheck() begin" << endl; - - ArithVar leaving = ARITHVAR_SENTINEL; - ArithVar entering = ARITHVAR_SENTINEL; - DeltaRational leavingShift = d_DELTA_ZERO; // The amount the leaving variable can change by without making the tableau inconsistent - DeltaRational leavingDelta = d_DELTA_ZERO; // The amount the optimization function changes by selecting leaving - - Assert(d_improvementCandidates.empty()); - - for( Tableau::RowIterator ri = d_tableau.basicRowIterator(d_optRow); !ri.atEnd(); ++ri){ - const Tableau::Entry& e = *ri; - - ArithVar curr = e.getColVar(); - if(curr == d_optRow){ continue; } - - - int sgn = e.getCoefficient().sgn(); - Assert(sgn != 0); - if( (sgn < 0 && d_partialModel.strictlyBelowUpperBound(curr)) || - (sgn > 0 && d_partialModel.strictlyAboveLowerBound(curr)) ){ - - d_improvementCandidates.push_back(&e); - } - } - - if(d_improvementCandidates.empty()){ - Debug("primal") << "primalCheck() end : global" << endl; - return GlobalMinimum; // No variable in the optimization function can be improved - } - - DeltaRational minimumShift; - DeltaRational currShift; - for(EntryVector::const_iterator ci = d_improvementCandidates.begin(), end_ci = d_improvementCandidates.end(); ci != end_ci; ++ci){ - const Tableau::Entry& e = *(*ci); - ArithVar curr = e.getColVar(); - - ArithVar currEntering; - bool progress; - - minimumShift = (leaving == ARITHVAR_SENTINEL ) ? leavingDelta/(e.getCoefficient().abs()) : d_DELTA_ZERO; - int sgn = e.getCoefficient().sgn(); - computeShift(curr, sgn < 0, progress, currEntering, currShift, minimumShift); - - if(currEntering == ARITHVAR_SENTINEL){ - d_improvementCandidates.clear(); - - Debug("primal") << "primalCheck() end : unbounded" << endl; - d_primalCarry.d_unbounded = curr; - return FoundUnboundedVariable; - }else if(progress) { - leaving = curr; - leavingShift = currShift; - leavingDelta = currShift * e.getCoefficient(); - entering = currEntering; - - Assert(leavingDelta < d_DELTA_ZERO); - - const int RECHECK_PERIOD = 10; - if(d_pivotsSinceErrorProgress % RECHECK_PERIOD != 0){ - // we can make progress, stop - break; - } - } - } - - if(leaving == ARITHVAR_SENTINEL){ - cout << "Nothing can make progress " << endl; - - const uint32_t THRESHOLD = 20; - if(d_pivotsSinceOptProgress <= THRESHOLD){ - - int index = rand() % d_improvementCandidates.size(); - leaving = (*d_improvementCandidates[index]).getColVar(); - entering = selectFirstValid(leaving, (*d_improvementCandidates[index]).getCoefficient().sgn() < 0); - }else{ // Bland's rule - bool increasing; - for(EntryVector::const_iterator ci = d_improvementCandidates.begin(), end_ci = d_improvementCandidates.end(); ci != end_ci; ++ci){ - const Tableau::Entry& e = *(*ci); - ArithVar curr = e.getColVar(); - leaving = (leaving == ARITHVAR_SENTINEL) ? curr : minVarOrder(*this, curr, leaving); - if(leaving == curr){ - increasing = (e.getCoefficient().sgn() < 0); - } - } - - entering = selectMinimumValid(leaving, increasing); - } - Assert(leaving != ARITHVAR_SENTINEL); - Assert(entering != ARITHVAR_SENTINEL); - - d_primalCarry.d_leaving = leaving; - d_primalCarry.d_entering = entering; - - d_primalCarry.d_nextEnteringValue = d_partialModel.getAssignment(entering); - - Debug("primal") << "primalCheck() end : no progress made " << leaving << " to " << entering << " (" << d_pivotsSinceOptProgress << ")"<< endl; - d_improvementCandidates.clear(); - return NoProgressOnLeaving; - }else{ - const Tableau::Entry& enterLeavingEntry = d_tableau.findEntry(d_tableau.basicToRowIndex(entering), leaving); - Assert(!enterLeavingEntry.blank()); - - d_primalCarry.d_leaving = leaving; - d_primalCarry.d_entering = entering; - d_primalCarry.d_nextEnteringValue = d_partialModel.getAssignment(entering) - + leavingShift * enterLeavingEntry.getCoefficient(); - - Debug("primal") << "primalCheck() end : progress" << endl - << leaving << " to " << entering << " ~ " - << d_partialModel.getAssignment(leaving) << " ~ " << leavingShift - << " ~ " << enterLeavingEntry.getCoefficient() - << " ~ " << d_primalCarry.d_nextEnteringValue << endl; - - d_improvementCandidates.clear(); - return MakeProgressOnLeaving; - } - - // anyProgress = true; - - // DeltaRational currDelta = currShift * e.getCoefficient(); - - // int cmp = currDelta.cmp(leavingDelta); - - // // Cases: - // // 1) No candidate yet, - // // 2) there is a candidate with a strictly better update, or - // // 3) there is a candidate with the same update value that has a smaller value in the variable ordering. - // // - // // Case 3 covers Bland's rule. - // if(entering == ARITHVAR_SENTINEL || cmp < 0){ - // leaving = curr; - // }else if( cmp == 0 ){ - // leaving = minVarOrder(*this, curr, leaving); - // } - - // if(leaving == curr){ - // leavingShift = currShift; - // leavingDelta = currDelta; - // entering = currEntering; - // } - // } - // } - - // if(leaving == ARITHVAR_SENTINEL){ - // Debug("primal") << "primalCheck() end : global" << endl; - // return GlobalMinimum; // No variable in the optimization function can be improved - // }else{ - // const Tableau::Entry& enterLeavingEntry = d_tableau.findEntry(d_tableau.basicToRowIndex(entering), leaving); - // Assert(!enterLeavingEntry.blank()); - - // d_primalCarry.d_leaving = leaving; - // d_primalCarry.d_entering = entering; - // d_primalCarry.d_nextEnteringValue = d_partialModel.getAssignment(entering) - // + leavingShift * enterLeavingEntry.getCoefficient(); - - // Debug("primal") << "primalCheck() end : progress" << endl - // << leaving << " to " << entering << " ~ " - // << d_partialModel.getAssignment(leaving) << " ~ " << leavingShift - // << " ~ " << enterLeavingEntry.getCoefficient() - // << " ~ " << d_primalCarry.d_nextEnteringValue << endl; - // return MakeProgressOnLeaving; - // } -} - -ArithVar SimplexDecisionProcedure::selectMinimumValid(ArithVar v, bool increasing){ - ArithVar minimum = ARITHVAR_SENTINEL; - for(Tableau::ColIterator colIter = d_tableau.colIterator(v);!colIter.atEnd(); ++colIter){ - const Tableau::Entry& e = *colIter; - ArithVar basic = d_tableau.rowIndexToBasic(e.getRowIndex()); - if(basic == d_optRow) continue; - - - int esgn = e.getCoefficient().sgn(); - bool basicInc = (increasing == (esgn > 0)); - - if(!(basicInc ? d_partialModel.strictlyBelowUpperBound(basic) : - d_partialModel.strictlyAboveLowerBound(basic))){ - if(minimum == ARITHVAR_SENTINEL){ - minimum = basic; - }else{ - minimum = minVarOrder(*this, basic, minimum); - } - } - } - return minimum; -} - -ArithVar SimplexDecisionProcedure::selectFirstValid(ArithVar v, bool increasing){ - ArithVar minimum = ARITHVAR_SENTINEL; - - for(Tableau::ColIterator colIter = d_tableau.colIterator(v);!colIter.atEnd(); ++colIter){ - const Tableau::Entry& e = *colIter; - ArithVar basic = d_tableau.rowIndexToBasic(e.getRowIndex()); - if(basic == d_optRow) continue; - - int esgn = e.getCoefficient().sgn(); - bool basicInc = (increasing == (esgn > 0)); - - if(!(basicInc ? d_partialModel.strictlyBelowUpperBound(basic) : - d_partialModel.strictlyAboveLowerBound(basic))){ - if(minimum == ARITHVAR_SENTINEL){ - minimum = basic; - }else{ - minimum = minRowLength(*this, basic, minimum); - } - } - } - return minimum; -} - - - -void SimplexDecisionProcedure::computeShift(ArithVar leaving, bool increasing, bool& progress, ArithVar& entering, DeltaRational& shift, const DeltaRational& minimumShift){ - Assert(increasing ? (minimumShift >= d_DELTA_ZERO) : (minimumShift <= d_DELTA_ZERO) ); - - static int instance = 0; - Debug("primal") << "computeshift " << ++instance << " " << leaving << endl; - - // The selection for the leaving variable - entering = ARITHVAR_SENTINEL; - - // no progress is initially made - progress = false; - - bool bounded = false; - - if(increasing ? d_partialModel.hasUpperBound(leaving) : d_partialModel.hasLowerBound(leaving)){ - const DeltaRational& assignment = d_partialModel.getAssignment(leaving); - - bounded = true; - - DeltaRational diff = increasing ? d_partialModel.getUpperBound(leaving) - assignment : d_partialModel.getLowerBound(leaving) - assignment; - Assert(increasing ? diff.sgn() >=0 : diff.sgn() <= 0); - if((increasing) ? (diff < minimumShift) : ( diff > minimumShift)){ - Assert(!progress); - entering = leaving; // My my my, what an ugly hack - return; // no progress is possible stop - } - } - - // shift has a meaningful value once entering has a meaningful value - // if increasing, - // then shift > minimumShift >= 0 - // else shift < minimumShift <= 0 - // - // Maintain the following invariant: - // - // if increasing, - // if e_ij > 0, diff >= shift > minimumShift >= 0 - // if e_ij < 0, diff >= shift > minimumShift >= 0 - // if !increasing, - // if e_ij > 0, diff <= shift < minimumShift <= 0 - // if e_ij < 0, diff <= shift < minimumShift <= 0 - // if increasing == (e_ij > 0), diff = (d_partialModel.getUpperBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient() - // if increasing != (e_ij > 0), diff = (d_partialModel.getLowerBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient() - - for(Tableau::ColIterator colIter = d_tableau.colIterator(leaving);!colIter.atEnd(); ++colIter){ - const Tableau::Entry& e = *colIter; - - ArithVar basic = d_tableau.rowIndexToBasic(e.getRowIndex()); - if(basic == d_optRow) continue; - - int esgn = e.getCoefficient().sgn(); - bool basicInc = (increasing == (esgn > 0)); - // If both are true, increasing the variable entering increases the basic variable - // If both are false, the entering variable is decreasing, but the coefficient is negative and the basic variable is increasing - // If exactly one is false, the basic variable is decreasing. - - Debug("primal::shift") << basic << " " << d_partialModel.hasUpperBound(basic) << " " - << d_partialModel.hasLowerBound(basic) << " " - << e.getCoefficient() << endl; - - if( (basicInc && d_partialModel.hasUpperBound(basic))|| - (!basicInc && d_partialModel.hasLowerBound(basic))){ - - if(!(basicInc ? d_partialModel.strictlyBelowUpperBound(basic) : - d_partialModel.strictlyAboveLowerBound(basic))){ - // diff == 0, as diff > minimumShift >= 0 or diff < minimumShift <= 0 - Assert(!progress); - entering = basic; - return; - }else{ - DeltaRational diff = basicInc ? - (d_partialModel.getUpperBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient() : - (d_partialModel.getLowerBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient(); - - if( entering == ARITHVAR_SENTINEL ){ - if(increasing ? (diff <= minimumShift) : (diff >= minimumShift)){ - Assert(!progress); - entering = basic; - return; - }else{ - Assert(increasing ? (diff > minimumShift) : (diff < minimumShift)); - shift = diff; - entering = basic; - bounded = true; - } - }else{ - if( increasing ? (diff < shift) : diff > shift){ - // a new min for increasing - // a new max for decreasing - - if(increasing ? (diff <= minimumShift) : (diff >= minimumShift)){ - Assert(!progress); - entering = basic; - return; - }else{ - Assert(increasing ? (diff > minimumShift) : (diff < minimumShift)); - shift = diff; - entering = basic; - } - } - } - } - } - } - - if(!bounded){ - // A totally unbounded variable - Assert(entering == ARITHVAR_SENTINEL); - progress = true; - return; - }else if(entering == ARITHVAR_SENTINEL){ - // We have a variable that is bounded only by its maximum - for(Tableau::ColIterator colIter = d_tableau.colIterator(leaving);!colIter.atEnd(); ++colIter){ - const Tableau::Entry& e = *colIter; - - ArithVar basic = d_tableau.rowIndexToBasic(e.getRowIndex()); - if(basic == d_optRow){ - continue; - } else{ - entering = basic; - break; - } - } - Assert(entering != ARITHVAR_SENTINEL); - - Assert(increasing ? d_partialModel.hasUpperBound(leaving) : d_partialModel.hasLowerBound(leaving)); - - const DeltaRational& assignment = d_partialModel.getAssignment(leaving); - DeltaRational diff = increasing ? d_partialModel.getUpperBound(leaving) - assignment : d_partialModel.getLowerBound(leaving) - assignment; - - shift = diff; - - Assert(increasing ? shift.sgn() >=0 : shift.sgn() <= 0); - Assert(increasing ? shift > minimumShift : shift < minimumShift); - - progress = true; - return; - }else{ - Assert(bounded); - progress = true; - - if((increasing ? d_partialModel.hasUpperBound(leaving) : d_partialModel.hasLowerBound(leaving) )){ - Assert(entering != ARITHVAR_SENTINEL); - const DeltaRational& assignment = d_partialModel.getAssignment(leaving); - DeltaRational diff = increasing ? d_partialModel.getUpperBound(leaving) - assignment : d_partialModel.getLowerBound(leaving) - assignment; - if((increasing) ? (diff < shift) : ( diff > shift)){ - shift = diff; - } - } - - Assert(increasing ? shift.sgn() >=0 : shift.sgn() <= 0); - Assert(increasing ? shift > minimumShift : shift < minimumShift); - return; - } - - - // if(! bounded || - // (increasing && diff < shift) || // a new min for increasing - // (!increasing && diff > shift)){ // a new max for decreasing - // bounded = true; - // shift = diff; - // entering = basic; - // } - // } - - // if(notAtTheBound && !blandMode){ - // DeltaRational diff = basicInc ? - // (d_partialModel.getUpperBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient() : - // (d_partialModel.getLowerBound(basic) - d_partialModel.getAssignment(basic)) / e.getCoefficient(); - - // if(! bounded || - // (increasing && diff < shift) || // a new min for increasing - // (!increasing && diff > shift)){ // a new max for decreasing - // bounded = true; - // shift = diff; - // entering = basic; - // } - // }else if (!notAtTheBound) { // basic is already exactly at the bound - // if(!blandMode){ // Enter into using Bland's rule - // blandMode = true; - // bounded = true; - // shift = d_DELTA_ZERO; - // entering = basic; - // }else{ - // entering = minVarOrder(*this, entering, basic); // Bland's rule. - // } - // } - // else this basic variable cannot be violated by increasing/decreasing entering - - - - - // if(!blandMode && (increasing ? d_partialModel.hasUpperBound(leaving) : d_partialModel.hasLowerBound(leaving) )){ - // Assert(entering != ARITHVAR_SENTINEL); - // bounded = true; - // DeltaRational diff = increasing ? d_partialModel.getUpperBound(leaving) - assignment : d_partialModel.getLowerBound(leaving) - assignment; - // if((increasing) ? (diff < shift) : ( diff > shift)){ - // shift = diff; - // } - // } - - // Assert(increasing ? shift.sgn() >=0 : shift.sgn() <= 0); - - // return shift; -} - - -/** - * Given an variable on the optimization row that can be used to decrease the value of the optimization function - * arbitrarily and an optimization function that is strictly positive in the current model, - * driveOptToZero updates the value of unbounded s.t. the value of d_opt is exactly 0. - */ -void SimplexDecisionProcedure::driveOptToZero(ArithVar unbounded){ - Assert(!belowThreshold()); - - const Tableau::Entry& e = d_tableau.findEntry(d_tableau.basicToRowIndex(d_optRow), unbounded); - Assert(!e.blank()); - - DeltaRational theta = (d_negOptConstant-d_partialModel.getAssignment(d_optRow))/ (e.getCoefficient()); - Assert((e.getCoefficient().sgn() > 0) ? (theta.sgn() < 0) : (theta.sgn() > 0)); - - DeltaRational newAssignment = d_partialModel.getAssignment(unbounded) + theta; - d_linEq.update(unbounded, newAssignment); - - if(Debug.isOn("primal::consistent")){ Assert(d_linEq.debugEntireLinEqIsConsistent("driveOptToZero")); } - - Assert(belowThreshold()); -} |