diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index aab35da90..bfd80e3af 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -331,7 +331,7 @@ pair> simplex(vector _constraints, Linear /// Turns all bounds into constraints. /// @returns false if the bounds make the state infeasible. -bool boundsToConstraints(SolvingState& _state) +optional boundsToConstraints(SolvingState& _state) { size_t columns = _state.variableNames.size(); @@ -341,14 +341,14 @@ bool boundsToConstraints(SolvingState& _state) if (bounds.lower && bounds.upper) { if (*bounds.lower > *bounds.upper) - return false; + return bounds.lowerReasons + bounds.upperReasons; if (*bounds.lower == *bounds.upper) { LinearExpression c; c.resize(columns); c[0] = *bounds.lower; c[varIndex] = bigint(1); - _state.constraints.emplace_back(Constraint{move(c), true}); + _state.constraints.emplace_back(Constraint{move(c), true, bounds.lowerReasons + bounds.upperReasons}); continue; } } @@ -358,7 +358,7 @@ bool boundsToConstraints(SolvingState& _state) c.resize(columns); c[0] = -*bounds.lower; c[varIndex] = bigint(-1); - _state.constraints.emplace_back(Constraint{move(c), false}); + _state.constraints.emplace_back(Constraint{move(c), false, move(bounds.lowerReasons)}); } if (bounds.upper) { @@ -366,11 +366,11 @@ bool boundsToConstraints(SolvingState& _state) c.resize(columns); c[0] = *bounds.upper; c[varIndex] = bigint(1); - _state.constraints.emplace_back(Constraint{move(c), false}); + _state.constraints.emplace_back(Constraint{move(c), false, move(bounds.upperReasons)}); } } _state.bounds.clear(); - return true; + return nullopt; } /// Removes incides set to true from a vector-like data structure. @@ -509,6 +509,7 @@ string SolvingState::toString() const { string result; + // TODO print reasons for (Constraint const& constraint: constraints) { vector line; @@ -544,27 +545,27 @@ string SolvingState::toString() const return result; } -pair> SolvingStateSimplifier::simplify() +pair> SolvingStateSimplifier::simplify() { do { m_changed = false; - if ( - !removeFixedVariables() || - !extractDirectConstraints() || - // Used twice on purpose - !removeFixedVariables() || - !removeEmptyColumns() - ) - return {LPResult::Infeasible, {}}; + if (auto conflict = removeFixedVariables()) + return {LPResult::Infeasible, move(*conflict)}; + if (auto conflict = extractDirectConstraints()) + return {LPResult::Infeasible, move(*conflict)}; + // Used twice on purpose + if (auto conflict = removeFixedVariables()) + return {LPResult::Infeasible, move(*conflict)}; + removeEmptyColumns(); } while (m_changed); return {LPResult::Unknown, move(m_model)}; } -bool SolvingStateSimplifier::removeFixedVariables() +optional SolvingStateSimplifier::removeFixedVariables() { for (auto const& [index, bounds]: m_state.bounds | ranges::views::enumerate) { @@ -574,26 +575,31 @@ bool SolvingStateSimplifier::removeFixedVariables() rational lower = max(rational{}, bounds.lower ? *bounds.lower : rational{}); rational upper = *bounds.upper; if (upper < lower) - return false; // Infeasible. + // Infeasible. + return bounds.lowerReasons + bounds.upperReasons; if (upper != lower) continue; + set reasons = bounds.lowerReasons + bounds.upperReasons; m_model[m_state.variableNames.at(index)] = lower; m_state.bounds[index] = {}; m_changed = true; // substitute variable + // -> add the bounds to the literals for the conflict + // (maybe only if one of these constraints is used) for (Constraint& constraint: m_state.constraints) if (constraint.data[index]) { constraint.data[0] -= constraint.data[index] * lower; constraint.data[index] = 0; + constraint.reasons += reasons; } } - return true; + return nullopt; } -bool SolvingStateSimplifier::extractDirectConstraints() +optional SolvingStateSimplifier::extractDirectConstraints() { vector constraintsToRemove(m_state.constraints.size(), false); bool needsRemoval = false; @@ -616,7 +622,7 @@ bool SolvingStateSimplifier::extractDirectConstraints() constraint.data.front().numerator() < 0 || (constraint.equality && constraint.data.front()) ) - return false; // Infeasible. + return constraint.reasons; } else { @@ -627,13 +633,19 @@ bool SolvingStateSimplifier::extractDirectConstraints() (factor >= 0 || constraint.equality) && (!m_state.bounds[varIndex].upper || bound < m_state.bounds[varIndex].upper) ) + { m_state.bounds[varIndex].upper = bound; + m_state.bounds[varIndex].upperReasons = constraint.reasons; + } if ( (factor <= 0 || constraint.equality) && + bound >= 0 && (!m_state.bounds[varIndex].lower || bound > m_state.bounds[varIndex].lower) ) - // Lower bound must be at least zero. - m_state.bounds[varIndex].lower = max(rational{}, bound); + { + m_state.bounds[varIndex].lower = bound; + m_state.bounds[varIndex].lowerReasons = constraint.reasons; + } } } if (needsRemoval) @@ -641,10 +653,10 @@ bool SolvingStateSimplifier::extractDirectConstraints() m_changed = true; eraseIndices(m_state.constraints, constraintsToRemove); } - return true; + return nullopt; } -bool SolvingStateSimplifier::removeEmptyColumns() +void SolvingStateSimplifier::removeEmptyColumns() { vector variablesSeen(m_state.bounds.size(), false); for (auto const& constraint: m_state.constraints) @@ -679,7 +691,6 @@ bool SolvingStateSimplifier::removeEmptyColumns() m_changed = true; removeColumns(m_state, variablesToRemove); } - return true; } SolvingState ProblemSplitter::next() @@ -722,7 +733,7 @@ SolvingState ProblemSplitter::next() // Use const& on purpose because moving from the state can lead // to undefined behaviour for connectedComponent Constraint const& constraint = m_state.constraints[i]; - Constraint splitRow{{}, constraint.equality}; + Constraint splitRow{{}, constraint.equality, constraint.reasons}; for (size_t j = 0; j < constraint.data.size(); j++) if (j == 0 || includedColumns[j]) splitRow.data.push_back(constraint.data[j]); @@ -738,15 +749,15 @@ LPSolver::LPSolver(bool _supportModels): { } -pair> LPSolver::check(SolvingState _state) +pair> LPSolver::check(SolvingState _state) { normalizeRowLengths(_state); - auto&& [simplificationResult, model] = SolvingStateSimplifier{_state}.simplify(); + auto&& [simplificationResult, modelOrReasonSet] = SolvingStateSimplifier{_state}.simplify(); switch (simplificationResult) { case LPResult::Infeasible: - return {LPResult::Infeasible, {}}; + return {LPResult::Infeasible, modelOrReasonSet}; case LPResult::Feasible: case LPResult::Unbounded: solAssert(false); @@ -754,6 +765,8 @@ pair> LPSolver::check(SolvingState _state) break; } + Model model = get(modelOrReasonSet); + bool canOnlyBeUnknown = false; ProblemSplitter splitter(move(_state)); while (splitter) @@ -765,25 +778,23 @@ pair> LPSolver::check(SolvingState _state) LPResult lpResult; vector solution; + if (auto conflict = boundsToConstraints(split)) + return {LPResult::Infeasible, move(*conflict)}; + auto it = m_cache.find(split); if (it != m_cache.end()) tie(lpResult, solution) = it->second; else { - SolvingState orig = split; - if (!boundsToConstraints(split)) - lpResult = LPResult::Infeasible; - else - { - LinearExpression objectives; - objectives.resize(1); - objectives.resize(split.constraints.front().data.size(), rational(bigint(1))); - tie(lpResult, solution) = simplex(split.constraints, move(objectives)); - } + LinearExpression objectives; + objectives.resize(1); + objectives.resize(split.constraints.front().data.size(), rational(bigint(1))); + tie(lpResult, solution) = simplex(split.constraints, move(objectives)); + // If we do not support models, do not store it in the cache because // the variable associations will be wrong. // Otherwise, it is fine to use the model. - m_cache.emplace(move(orig), make_pair(lpResult, m_supportModels ? solution : vector{})); + m_cache.emplace(split, make_pair(lpResult, m_supportModels ? solution : vector{})); } switch (lpResult) @@ -792,7 +803,13 @@ pair> LPSolver::check(SolvingState _state) case LPResult::Unbounded: break; case LPResult::Infeasible: - return {LPResult::Infeasible, {}}; + { + solAssert(split.bounds.empty()); + set reasons; + for (auto const& constraint: split.constraints) + reasons += constraint.reasons; + return {LPResult::Infeasible, move(reasons)}; + } case LPResult::Unknown: // We do not stop here, because another independent query can still be infeasible. canOnlyBeUnknown = true; @@ -804,7 +821,7 @@ pair> LPSolver::check(SolvingState _state) } if (canOnlyBeUnknown) - return {LPResult::Unknown, {}}; + return {LPResult::Unknown, Model{}}; return {LPResult::Feasible, move(model)}; } diff --git a/libsolutil/LP.h b/libsolutil/LP.h index 9fd50aca1..e5c11c2ec 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -28,6 +28,9 @@ namespace solidity::util { +using Model = std::map; +using ReasonSet = std::set; + /** * Constraint of the form * - data[1] * x_1 + data[2] * x_2 + ... <= data[0] (equality == false) @@ -38,6 +41,8 @@ struct Constraint { LinearExpression data; bool equality = false; + /// Set of literals the conjunction of which implies this constraint. + std::set reasons; bool operator<(Constraint const& _other) const; bool operator==(Constraint const& _other) const; @@ -57,10 +62,17 @@ struct SolvingState std::optional upper; bool operator<(Bounds const& _other) const { return make_pair(lower, upper) < make_pair(_other.lower, _other.upper); } bool operator==(Bounds const& _other) const { return make_pair(lower, upper) == make_pair(_other.lower, _other.upper); } + + /// Set of literals the conjunction of which implies the lower bonud. + std::set lowerReasons; + /// Set of literals the conjunction of which implies the upper bonud. + std::set upperReasons; }; /// Lower and upper bounds for variables (in the sense of >= / <=). std::vector bounds; std::vector constraints; + // For each bound and constraint, store an index of the literal + // that implies it. struct Compare { @@ -103,25 +115,26 @@ public: SolvingStateSimplifier(SolvingState& _state): m_state(_state) {} - std::pair> simplify(); + std::pair> simplify(); private: /// Remove variables that have equal lower and upper bound. - /// @returns false if the system is infeasible. - bool removeFixedVariables(); + /// @returns reason / set of conflicting clauses if infeasible. + std::optional removeFixedVariables(); /// Removes constraints of the form 0 <= b or 0 == b (no variables) and /// turns constraints of the form a * x <= b (one variable) into bounds. - bool extractDirectConstraints(); + /// @returns reason / set of conflicting clauses if infeasible. + std::optional extractDirectConstraints(); /// Removes all-zeros columns. - bool removeEmptyColumns(); + void removeEmptyColumns(); /// Set to true by the strategies if they performed some changes. bool m_changed = false; SolvingState& m_state; - std::map m_model; + Model m_model; }; /** @@ -167,7 +180,7 @@ class LPSolver public: explicit LPSolver(bool _supportModels = true); - std::pair>> check(SolvingState _state); + std::pair> check(SolvingState _state); private: using CacheValue = std::pair>>; diff --git a/test/libsolutil/LP.cpp b/test/libsolutil/LP.cpp index 8bc6edd57..c88b0735d 100644 --- a/test/libsolutil/LP.cpp +++ b/test/libsolutil/LP.cpp @@ -52,11 +52,11 @@ public: } /// Adds the constraint "_lhs <= _rhs". - void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs) + void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs, set _reason = {}) { _lhs -= _rhs; _lhs[0] = -_lhs[0]; - m_solvingState.constraints.push_back({move(_lhs), false}); + m_solvingState.constraints.push_back({move(_lhs), false, move(_reason)}); } void addLEConstraint(LinearExpression _lhs, rational _rhs) @@ -65,33 +65,36 @@ public: } /// Adds the constraint "_lhs = _rhs". - void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs) + void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs, set _reason = {}) { _lhs -= _rhs; _lhs[0] = -_lhs[0]; - m_solvingState.constraints.push_back({move(_lhs), true}); + m_solvingState.constraints.push_back({move(_lhs), true, move(_reason)}); } - void addLowerBound(string _variable, rational _value) + void addLowerBound(string _variable, rational _value, set _reason = {}) { size_t index = variableIndex(_variable); if (index >= m_solvingState.bounds.size()) m_solvingState.bounds.resize(index + 1); m_solvingState.bounds.at(index).lower = _value; + m_solvingState.bounds.at(index).lowerReasons = move(_reason); } - void addUpperBound(string _variable, rational _value) + void addUpperBound(string _variable, rational _value, set _reason = {}) { size_t index = variableIndex(_variable); if (index >= m_solvingState.bounds.size()) m_solvingState.bounds.resize(index + 1); m_solvingState.bounds.at(index).upper = _value; + m_solvingState.bounds.at(index).upperReasons = move(_reason); } void feasible(vector> const& _solution) { - auto [result, model] = m_solver.check(m_solvingState); + auto [result, modelOrReasonSet] = m_solver.check(m_solvingState); BOOST_REQUIRE(result == LPResult::Feasible); + Model model = get(modelOrReasonSet); for (auto const& [var, value]: _solution) BOOST_CHECK_MESSAGE( value == model.at(var), @@ -102,6 +105,7 @@ public: void infeasible() { auto [result, model] = m_solver.check(m_solvingState); + // TODO check reason set BOOST_CHECK(result == LPResult::Infeasible); } @@ -283,7 +287,7 @@ BOOST_AUTO_TEST_CASE(equal_constant) feasible({{"x", 5}, {"y", 5}}); } -BOOST_AUTO_TEST_CASE(equal_constant) +BOOST_AUTO_TEST_CASE(all_equality) { auto x = variable("x"); auto y = variable("y");