From 2ef255bdab1b213b6208dd0943c55386262112f7 Mon Sep 17 00:00:00 2001 From: chriseth Date: Fri, 1 Apr 2022 19:01:50 +0200 Subject: [PATCH] Turn into class. --- libsolutil/LP.cpp | 416 ++++++++++++++++++++++------------------------ libsolutil/LP.h | 24 +++ 2 files changed, 219 insertions(+), 221 deletions(-) diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index 1db7a3b0c..1dc17d553 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -381,225 +381,6 @@ pair> simplex(vector _constraints, Linear return make_pair(result, solutionVector(tableau)); } -/// Introduces a new variable for each now and returns -/// a map from the indices of those variables -/// to the constraint they correspond to. -map normalizeForSolvingWithBounds(SolvingState& _state) -{ - size_t varsNeeded = _state.constraints.size(); - /* - for (Constraint& c: _state.constraints) - if (!c.equality || c.data[0]) - varsNeeded++; - */ - - map basicVariables; - size_t row = 0; - for (Constraint& c: _state.constraints) - { - c.data.resize(_state.variableNames.size() + varsNeeded); - /* - if (c.equality && !c.data[0]) - continue; - */ - - // ax + by <= c - // -> ax + by - s = 0, 0 <= s <= c - // ax + by = c - // -> ax + by - s = 0, c <= s <= c - size_t newVarIndex = _state.variableNames.size(); - basicVariables[newVarIndex] = row++; - solAssert(_state.variableNames.size() == _state.bounds.size()); - // TODO name needed unique? - _state.variableNames.emplace_back("_s" + to_string(newVarIndex)); - _state.bounds.emplace_back(SolvingState::Bounds{ - {c.equality ? rational{c.data[0]} : rational{0}}, - {c.data[0]}, - {}, - {} - }); - c.equality = true; - solAssert(c.data.size() > newVarIndex); - c.data[newVarIndex] = -1; - c.data[0] = 0; - } - - return basicVariables; -} - -void withBoundsUpdate( - SolvingState& _state, - vector& _assignments, - map const& _basicVariables, - size_t _i, - rational const& _value -) -{ - rational delta = _value - _assignments[_i]; - _assignments[_i] = _value; - for (size_t j = 0; j < _assignments.size(); j++) - if (_basicVariables.count(j) && _state.constraints[_basicVariables.at(j)].data[_i]) - _assignments[j] += _state.constraints[_basicVariables.at(j)].data[_i] * delta; -} - -optional firstConflictingBasicVariable( - SolvingState& _state, - vector& _assignments, - map const& _basicVariables -) -{ - for (auto const& [i, bounds]: _state.bounds | ranges::views::enumerate) - if (_basicVariables.count(i) && ( - (bounds.lower && _assignments[i] < *bounds.lower) || - (bounds.upper && _assignments[i] > *bounds.upper) - )) - return i; - return nullopt; -} - -optional firstReplacementVar( - SolvingState& _state, - vector& _assignments, - map const& _basicVariables, - size_t _basicVarToReplace, - bool _increasing -) -{ - LinearExpression const& basicVarEquation = _state.constraints[_basicVariables.at(_basicVarToReplace)].data; - for (auto const& [i, bounds]: _state.bounds | ranges::views::enumerate) - { - if (_basicVariables.count(i) || !basicVarEquation[i]) - continue; - bool positive = basicVarEquation[i] > 0; - if (!_increasing) - positive = !positive; - if (positive && (!_state.bounds[i].upper || _assignments[i] < _state.bounds[i].upper)) - return i; - if (!positive && (!_state.bounds[i].lower || _assignments[i] > _state.bounds[i].lower)) - return i; - } - return nullopt; -} - -void pivot(SolvingState& _state, map& _basicVariables, size_t _old, size_t _new) -{ - // Transform pivotRow such that the coefficient for _new is -1 - // Then use that to set all other coefficients for _new to zero. - size_t pivotRow = _basicVariables[_old]; - LinearExpression& pivotRowData = _state.constraints[_basicVariables[_old]].data; - - rational pivot = pivotRowData[_new]; - solAssert(pivot != 0, ""); - if (pivot != -1) - pivotRowData /= -pivot; - solAssert(pivotRowData[_new] == rational(-1), ""); - - auto subtractMultipleOfPivotRow = [&](LinearExpression& _row) { - if (_row[_new] == 0) - return; - else if (_row[_new] == rational{1}) - _row += pivotRowData; - else if (_row[_new] == rational{-1}) - _row -= pivotRowData; - else - _row += _row[_new] * pivotRowData; - }; - - for (size_t i = 0; i < _state.constraints.size(); ++i) - if (i != pivotRow) - subtractMultipleOfPivotRow(_state.constraints[i].data); - - _basicVariables.erase(_old); - _basicVariables[_new] = pivotRow; -} - -void pivotAndUpdate( - SolvingState& _state, - map& _basicVariables, - vector& _assignments, - size_t _oldBasicVar, - rational const& _newValue, - size_t _newBasicVar -) -{ - rational theta = (_newValue - _assignments[_oldBasicVar]) / _state.constraints[_basicVariables[_oldBasicVar]].data[_newBasicVar]; - - _assignments[_oldBasicVar] = _newValue; - _assignments[_newBasicVar] += theta; - - for (auto const& [i, row]: _basicVariables) - if (i != _oldBasicVar && _state.constraints[row].data[_newBasicVar]) - _assignments[i] += _state.constraints[row].data[_newBasicVar] * theta; - - pivot(_state, _basicVariables, _oldBasicVar, _newBasicVar); - cout << "After pivot and update: " << endl << _state.toString() << endl; -} - - -pair> simplexWithBounds(SolvingState _state) -{ - cout << "===========================" << endl; - cout << "===========================" << endl; - cout << "simpl with bounds on\n" << _state.toString() << endl; - map basicVariables = normalizeForSolvingWithBounds(_state); - cout << "After norm ------------------------\n" << _state.toString() << endl; - - vector assignments(_state.variableNames.size()); - // We start with an all-zero assignment and then gradually add - // the bounds we already have. - - cout << "Adjusting bounds on non-basic variables\n"; - // Adjust the assignments so we satisfy the bounds of the non-basic variables. - for (auto const& [i, bounds]: _state.bounds | ranges::views::enumerate) - { - if (basicVariables.count(i) || (!bounds.lower && !bounds.upper)) - continue; - if (bounds.lower && bounds.upper) - solAssert(*bounds.lower <= *bounds.upper); - if (bounds.lower && assignments[i] < *bounds.lower) - withBoundsUpdate(_state, assignments, basicVariables, i, *bounds.lower); - else if (bounds.upper && assignments[i] > *bounds.upper) - withBoundsUpdate(_state, assignments, basicVariables, i, *bounds.upper); - cout << "Assignments after satisfying bound for " << _state.variableNames[i] << ":\n"; - for (auto const& [j, val]: assignments | ranges::views::enumerate) - { - cout << " - " << _state.variableNames[j]; - if (basicVariables.count(j)) cout << " (b)"; - cout << " = " << ::toString(val) << endl; - } - cout << "----------------\n" << _state.toString() << endl; - } - cout << "Bounds on non-basic vaiables set.\n"; - - // Now try to make the basic variables happy, pivoting if necessary. - - // TODO bound number of iterations - while (auto bvi = firstConflictingBasicVariable(_state, assignments, basicVariables)) - { - cout << "Basic variable " << _state.variableNames[*bvi] << " is conflicting." << endl; - if (_state.bounds[*bvi].lower && assignments[*bvi] < *_state.bounds[*bvi].lower) - { - if (auto replacementVar = firstReplacementVar(_state, assignments, basicVariables, *bvi, true)) - pivotAndUpdate(_state, basicVariables, assignments, *bvi, *_state.bounds.at(*bvi).lower, *replacementVar); - else - return make_pair(LPResult::Infeasible, vector{}); - } - else if (_state.bounds[*bvi].upper && assignments[*bvi] > *_state.bounds[*bvi].upper) - { - if (auto replacementVar = firstReplacementVar(_state, assignments, basicVariables, *bvi, false)) - pivotAndUpdate(_state, basicVariables, assignments, *bvi, *_state.bounds.at(*bvi).upper, *replacementVar); - else - return make_pair(LPResult::Infeasible, vector{}); - } - } - - cout << ">>>>>>>>>>>>>>>>>>>>>>" << endl; - cout << ">>>>>>>>>>>>>>>>>>>>>>" << endl; - - return make_pair(LPResult::Feasible, move(assignments)); -} - - /// Turns all bounds into constraints. /// @returns false if the bounds make the state infeasible. optional boundsToConstraints(SolvingState& _state) @@ -819,6 +600,199 @@ string SolvingState::toString() const return result; } +SimplexWithBounds::SimplexWithBounds(SolvingState _state): + m_state(move(_state)) +{ + size_t varsNeeded = m_state.constraints.size(); + size_t row = 0; + for (Constraint& c: m_state.constraints) + { + c.data.resize(m_state.variableNames.size() + varsNeeded); + size_t newVarIndex = m_state.variableNames.size(); + m_basicVariables[newVarIndex] = row++; + solAssert(m_state.variableNames.size() == m_state.bounds.size()); + // TODO name needed unique? + m_state.variableNames.emplace_back("_s" + to_string(newVarIndex)); + // TODO we assume zero as lower bound here, but we do not have to. + m_state.bounds.emplace_back(SolvingState::Bounds{ + {c.equality ? rational{c.data[0]} : rational{0}}, + {c.data[0]}, + {}, + {} + }); + c.equality = true; + solAssert(c.data.size() > newVarIndex); + c.data[newVarIndex] = -1; + c.data[0] = 0; + } + + m_assignments.resize(m_state.variableNames.size()); +} + +LPResult SimplexWithBounds::check() +{ +// cout << "===========================" << endl; +// cout << "===========================" << endl; +// cout << "simpl with bounds on\n" << toString() << endl; + + // We start with an all-zero assignment and then gradually add + // the bounds we already have. + + //cout << "Adjusting bounds on non-basic variables\n"; + // Adjust the assignments so we satisfy the bounds of the non-basic variables. + for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate) + { + if (m_basicVariables.count(i) || (!bounds.lower && !bounds.upper)) + continue; + if (bounds.lower && bounds.upper) + solAssert(*bounds.lower <= *bounds.upper); + if (bounds.lower && m_assignments[i] < *bounds.lower) + update(i, *bounds.lower); + else if (bounds.upper && m_assignments[i] > *bounds.upper) + update(i, *bounds.upper); + } + //cout << "Bounds on non-basic vaiables set.\n" << toString() << endl; + + // Now try to make the basic variables happy, pivoting if necessary. + + // TODO bound number of iterations + while (auto bvi = firstConflictingBasicVariable()) + { + //cout << "Basic variable " << m_state.variableNames[*bvi] << " is conflicting." << endl; + if (m_state.bounds[*bvi].lower && m_assignments[*bvi] < *m_state.bounds[*bvi].lower) + { + if (auto replacementVar = firstReplacementVar(*bvi, true)) + pivotAndUpdate(*bvi, *m_state.bounds.at(*bvi).lower, *replacementVar); + else + { + //cout << "------->>>>>>>>>>>>>> infeasible\n"; + return LPResult::Infeasible; + } + } + else if (m_state.bounds[*bvi].upper && m_assignments[*bvi] > *m_state.bounds[*bvi].upper) + { + if (auto replacementVar = firstReplacementVar(*bvi, false)) + pivotAndUpdate(*bvi, *m_state.bounds.at(*bvi).upper, *replacementVar); + else + { + //cout << "------->>>>>>>>>>>>>> infeasible\n"; + return LPResult::Infeasible; + } + } + } + + //cout << "SAT ------->>>>>>>>>>>>>>\n"; + + return LPResult::Feasible; + +} + +void SimplexWithBounds::update(size_t _var, rational const& _value) +{ + rational delta = _value - m_assignments[_var]; + m_assignments[_var] = _value; + for (size_t j = 0; j < m_assignments.size(); j++) + if (m_basicVariables.count(j) && m_state.constraints[m_basicVariables.at(j)].data[_var]) + m_assignments[j] += m_state.constraints[m_basicVariables.at(j)].data[_var] * delta; + +} + +optional SimplexWithBounds::firstConflictingBasicVariable() const +{ + for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate) + if (m_basicVariables.count(i) && ( + (bounds.lower && m_assignments[i] < *bounds.lower) || + (bounds.upper && m_assignments[i] > *bounds.upper) + )) + return i; + return nullopt; +} + +optional SimplexWithBounds::firstReplacementVar( + size_t _basicVarToReplace, + bool _increasing +) const +{ + LinearExpression const& basicVarEquation = m_state.constraints[m_basicVariables.at(_basicVarToReplace)].data; + for (auto const& [i, bounds]: m_state.bounds | ranges::views::enumerate) + { + if (m_basicVariables.count(i) || !basicVarEquation[i]) + continue; + bool positive = basicVarEquation[i] > 0; + if (!_increasing) + positive = !positive; + if (positive && (!m_state.bounds[i].upper || m_assignments[i] < m_state.bounds[i].upper)) + return i; + if (!positive && (!m_state.bounds[i].lower || m_assignments[i] > m_state.bounds[i].lower)) + return i; + } + return nullopt; +} + +void SimplexWithBounds::pivot(size_t _old, size_t _new) +{ + // Transform pivotRow such that the coefficient for _new is -1 + // Then use that to set all other coefficients for _new to zero. + size_t pivotRow = m_basicVariables[_old]; + LinearExpression& pivotRowData = m_state.constraints[m_basicVariables[_old]].data; + + rational pivot = pivotRowData[_new]; + solAssert(pivot != 0, ""); + if (pivot != -1) + pivotRowData /= -pivot; + solAssert(pivotRowData[_new] == rational(-1), ""); + + auto subtractMultipleOfPivotRow = [&](LinearExpression& _row) { + if (_row[_new] == 0) + return; + else if (_row[_new] == rational{1}) + _row += pivotRowData; + else if (_row[_new] == rational{-1}) + _row -= pivotRowData; + else + _row += _row[_new] * pivotRowData; + }; + + for (size_t i = 0; i < m_state.constraints.size(); ++i) + if (i != pivotRow) + subtractMultipleOfPivotRow(m_state.constraints[i].data); + + m_basicVariables.erase(_old); + m_basicVariables[_new] = pivotRow; +} + +void SimplexWithBounds::pivotAndUpdate( + size_t _oldBasicVar, + rational const& _newValue, + size_t _newBasicVar +) +{ + rational theta = (_newValue - m_assignments[_oldBasicVar]) / m_state.constraints[m_basicVariables[_oldBasicVar]].data[_newBasicVar]; + + m_assignments[_oldBasicVar] = _newValue; + m_assignments[_newBasicVar] += theta; + + for (auto const& [i, row]: m_basicVariables) + if (i != _oldBasicVar && m_state.constraints[row].data[_newBasicVar]) + m_assignments[i] += m_state.constraints[row].data[_newBasicVar] * theta; + + pivot(_oldBasicVar, _newBasicVar); + //cout << "After pivot and update: " << endl << toString() << endl; +} + +string SimplexWithBounds::toString() const +{ + string ret = "---\n" + m_state.toString(); + for (auto const& [j, val]: m_assignments | ranges::views::enumerate) + { + ret += " " + m_state.variableNames[j]; + if (m_basicVariables.count(j)) ret += " (b)"; + ret += " = " + ::toString(val) + "\n"; + } + ret += "------\n"; + return ret; +} + pair, ReasonSet>> SolvingStateSimplifier::simplify() { do @@ -1150,8 +1124,7 @@ pair> LPSolver::check() } if (!result) { - result = LPResult::Unknown; - tie(*result, problem->model) = simplexWithBounds(state); + *result = SimplexWithBounds{state}.check(); // TODO we should even keep the updated tableau in the subproblem if (m_cache) { @@ -1243,6 +1216,7 @@ void LPSolver::addConstraintToSubProblem(size_t _subProblem, Constraint _constra m_subProblemsPerVariable[index] = _subProblem; problem.variables.emplace(index); } + // TODO if we keep the matrix, we could just add this as a new row. problem.removableConstraints.emplace_back(move(_constraint)); } diff --git a/libsolutil/LP.h b/libsolutil/LP.h index 8e544e544..f750d6d52 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -172,6 +172,30 @@ enum class LPResult }; +class SimplexWithBounds +{ +public: + explicit SimplexWithBounds(SolvingState _state); + LPResult check(); + + std::string toString() const; +private: + /// Set value of non-basic variable. + void update(size_t _var, rational const& _value); + /// @returns the index of the first basic variable violating its bounds. + std::optional firstConflictingBasicVariable() const; + std::optional firstReplacementVar(size_t _basicVarToReplace, bool _increasing) const; + + void pivot(size_t _old, size_t _new); + void pivotAndUpdate(size_t _oldBasicVar, rational const& _newValue, size_t _newBasicVar); + + SolvingState m_state; + std::vector m_assignments; + /// Variable index to row it controls. + std::map m_basicVariables; +}; + + /** * Applies several strategies to simplify a given solving state. * During these simplifications, it can sometimes already be determined if the