diff --git a/libsolutil/BooleanLP.cpp b/libsolutil/BooleanLP.cpp index a5db6f658..9022a94fa 100644 --- a/libsolutil/BooleanLP.cpp +++ b/libsolutil/BooleanLP.cpp @@ -265,13 +265,13 @@ pair> BooleanLPSolver::check(vector cons Constraint const& constraint = state().conditionalConstraints.at(constraintIndex); lpSolvers.back().second.addConstraint(constraint, constraintIndex); } - auto&& [result, modelOrReason] = lpSolvers.back().second.check(); + auto&& [result, reasonSet] = lpSolvers.back().second.check(); // We can only really use the result "infeasible". Everything else should be "sat". if (result == LPResult::Infeasible) { // TODO is it ok to ignore the non-constraint boolean variables here? Clause conflictClause; - for (size_t constraintIndex: get(modelOrReason)) + for (size_t constraintIndex: reasonSet) conflictClause.emplace_back(Literal{false, constraintIndex}); return conflictClause; } diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index 5173dffd0..26e916b09 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -218,20 +218,28 @@ void LPSolver::setVariableName(size_t _variable, string _name) void LPSolver::addLowerBound(size_t _variable, rational _bound) { SubProblem& p = unsealForVariable(_variable); - Variable& var = p.variables[p.varMapping.at(_variable)]; + size_t innerIndex = p.varMapping.at(_variable); + Variable& var = p.variables[innerIndex]; if (!var.bounds.lower || *var.bounds.lower < _bound) + { var.bounds.lower = move(_bound); + p.variablesPotentiallyOutOfBounds.insert(innerIndex); + } } void LPSolver::addUpperBound(size_t _variable, rational _bound) { SubProblem& p = unsealForVariable(_variable); - Variable& var = p.variables[p.varMapping.at(_variable)]; + size_t innerIndex = p.varMapping.at(_variable); + Variable& var = p.variables[innerIndex]; if (!var.bounds.upper || *var.bounds.upper > _bound) + { var.bounds.upper = move(_bound); + p.variablesPotentiallyOutOfBounds.insert(innerIndex); + } } -pair> LPSolver::check() +pair LPSolver::check() { for (auto&& [index, problem]: m_subProblems | ranges::views::enumerate) if (problem) @@ -248,7 +256,7 @@ pair> LPSolver::check() return {LPResult::Infeasible, problem->reasons}; } //cout << "Feasible:\n" << toString() << endl; - return {LPResult::Feasible, model()}; + return {LPResult::Feasible, {}}; } string LPSolver::toString() const @@ -260,6 +268,16 @@ string LPSolver::toString() const return result; } +map LPSolver::model() const +{ + map result; + for (auto const& problem: m_subProblems) + if (problem) + for (auto&& [outerIndex, innerIndex]: problem->varMapping) + result[problem->variables[innerIndex].name] = problem->variables[innerIndex].value; + return result; +} + LPSolver::SubProblem& LPSolver::unseal(size_t _problemIndex) { shared_ptr& problem = m_subProblems[_problemIndex]; @@ -302,6 +320,8 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom) combineInto.factors.emplace_back(move(shiftedRow)); } combineInto.variables += combineFrom.variables; + for (auto const& index: combineFrom.variablesPotentiallyOutOfBounds) + combineInto.variablesPotentiallyOutOfBounds.insert(index + varShift); for (auto&& [index, row]: combineFrom.basicVariables) combineInto.basicVariables.emplace(index + varShift, row + rowShift); for (auto&& [outerIndex, innerIndex]: combineFrom.varMapping) @@ -372,7 +392,8 @@ void LPSolver::addConstraintToSubProblem( if (_constraint.equality) problem.variables[slackIndex].bounds.lower = _constraint.data[0]; problem.variables[slackIndex].bounds.upper = _constraint.data[0]; - + // TODO it is a basic var, so we don't add it, unless we use this for basic vars. + //problem.variablesPotentiallyOutOfBounds.insert(slackIndex); // Compress the constraint, i.e. turn outer variable indices into // inner variable indices. @@ -419,20 +440,6 @@ size_t LPSolver::addNewVariableToSubProblem(size_t _subProblem) return index; } -map LPSolver::model() const -{ - map result; - /* - * // This is actually unused - for (auto const& problem: m_subProblems) - if (problem) - for (auto&& [outerIndex, innerIndex]: problem->varMapping) - result[problem->variables[innerIndex].name] = problem->variables[innerIndex].value; - cout << "MOdel size: " << to_string(result.size()) << endl; - */ - return result; -} - LPResult LPSolver::SubProblem::check() { @@ -445,17 +452,8 @@ LPResult LPSolver::SubProblem::check() // cout << "----------------------------" << endl; // cout << "fixing non-basic..." << endl; // Adjust the assignments so we satisfy the bounds of the non-basic variables. - for (auto const& [i, var]: variables | ranges::views::enumerate) - { - if (var.bounds.lower && var.bounds.upper && *var.bounds.lower > *var.bounds.upper) - return LPResult::Infeasible; - if (basicVariables.count(i) || (!var.bounds.lower && !var.bounds.upper)) - continue; - if (var.bounds.lower && var.value < *var.bounds.lower) - update(i, *var.bounds.lower); - else if (var.bounds.upper && var.value > *var.bounds.upper) - update(i, *var.bounds.upper); - } + if (!correctNonbasic()) + return LPResult::Infeasible; // Now try to make the basic variables happy, pivoting if necessary. @@ -552,26 +550,54 @@ string LPSolver::SubProblem::toString() const return resultString + "----\n"; } +bool LPSolver::SubProblem::correctNonbasic() +{ + set toCorrect; + swap(toCorrect, variablesPotentiallyOutOfBounds); + for (size_t i: toCorrect) + { + Variable& var = variables.at(i); + if (var.bounds.lower && var.bounds.upper && *var.bounds.lower > *var.bounds.upper) + return false; + if (basicVariables.count(i)) + { + variablesPotentiallyOutOfBounds.insert(i); + continue; + } + if (!var.bounds.lower && !var.bounds.upper) + continue; + if (var.bounds.lower && var.value < *var.bounds.lower) + update(i, *var.bounds.lower); + else if (var.bounds.upper && var.value > *var.bounds.upper) + update(i, *var.bounds.upper); + } + return true; +} + void LPSolver::SubProblem::update(size_t _varIndex, rational const& _value) { rational delta = _value - variables[_varIndex].value; variables[_varIndex].value = _value; for (size_t j = 0; j < variables.size(); j++) if (basicVariables.count(j) && factors[basicVariables.at(j)][_varIndex]) + { variables[j].value += delta * factors[basicVariables.at(j)][_varIndex]; + //variablesPotentiallyOutOfBounds.insert(j); + } } optional LPSolver::SubProblem::firstConflictingBasicVariable() const { - for (auto const& varItem: basicVariables) + // TODO we could use "variablesPotentiallyOutOfBounds" here. + for (auto&& [i, row]: basicVariables) { - Variable const& variable = variables[varItem.first]; + Variable const& variable = variables[i]; if ( (variable.bounds.lower && variable.value < *variable.bounds.lower) || (variable.bounds.upper && variable.value > *variable.bounds.upper) ) - return varItem.first; + return i; } return nullopt; } @@ -582,14 +608,14 @@ optional LPSolver::SubProblem::firstReplacementVar( ) const { LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)]; - for (auto const& [i, var]: variables | ranges::views::enumerate) + for (auto const& [i, factor]: basicVarEquation.enumerate()) { - if (basicVariables.count(i) || !basicVarEquation[i]) + if (i == _basicVarToReplace || !factor) continue; - bool positive = basicVarEquation[i] > 0; + bool positive = factor > 0; if (!_increasing) positive = !positive; - Variable const& candidate = variables[i]; + Variable const& candidate = variables.at(i); if (positive && (!candidate.bounds.upper || candidate.value < *candidate.bounds.upper)) return i; if (!positive && (!candidate.bounds.lower || candidate.value > *candidate.bounds.lower)) diff --git a/libsolutil/LP.h b/libsolutil/LP.h index c1ea8734d..0b4e9693f 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -196,9 +196,10 @@ public: void addLowerBound(size_t _variable, rational _bound); void addUpperBound(size_t _variable, rational _bound); - std::pair> check(); + std::paircheck(); std::string toString() const; + std::map model() const; private: struct Bounds @@ -219,6 +220,7 @@ private: std::optional result = std::nullopt; std::vector factors; std::vector variables; + std::set variablesPotentiallyOutOfBounds; /// Variable index to constraint it controls. std::map basicVariables; /// Maps outer indices to inner indices. @@ -228,6 +230,7 @@ private: LPResult check(); std::string toString() const; private: + bool correctNonbasic(); /// Set value of non-basic variable. void update(size_t _varIndex, rational const& _value); /// @returns the index of the first basic variable violating its bounds. @@ -247,8 +250,6 @@ private: void addOuterVariableToSubProblem(size_t _subProblem, size_t _outerIndex); size_t addNewVariableToSubProblem(size_t _subProblem); - std::map model() const; - /// These use "copy on write". std::vector> m_subProblems; /// Maps outer indices to sub problems. diff --git a/test/libsolutil/LP.cpp b/test/libsolutil/LP.cpp index 331190479..efc9956b9 100644 --- a/test/libsolutil/LP.cpp +++ b/test/libsolutil/LP.cpp @@ -83,9 +83,9 @@ public: void feasible(vector> const& _solution) { - auto [result, modelOrReasonSet] = m_solver.check(); + auto [result, reasonSet] = m_solver.check(); BOOST_REQUIRE(result == LPResult::Feasible); - Model model = get(modelOrReasonSet); + Model model = m_solver.model(); for (auto const& [var, value]: _solution) { BOOST_CHECK_MESSAGE( @@ -101,10 +101,9 @@ public: void infeasible(set _reason = {}) { - auto [result, modelOrReason] = m_solver.check(); + auto [result, reasonSet] = m_solver.check(); BOOST_REQUIRE(result == LPResult::Infeasible); - ReasonSet suppliedReason = get(modelOrReason); - BOOST_CHECK_MESSAGE(suppliedReason == _reason, "Reasons are different"); + BOOST_CHECK_MESSAGE(reasonSet == _reason, "Reasons are different"); } protected: