From 5faebbff39a81d42ce6a5e348cfde9a0aaf236b4 Mon Sep 17 00:00:00 2001 From: chriseth Date: Thu, 3 Feb 2022 13:12:50 +0100 Subject: [PATCH] Extract simplification class. --- libsolutil/LP.cpp | 310 ++++++++++++++++++++++------------------------ libsolutil/LP.h | 50 +++++++- 2 files changed, 197 insertions(+), 163 deletions(-) diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index a1d6ae8a2..53103c886 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -394,133 +394,12 @@ void removeColumns(SolvingState& _state, vector const& _columnsToRemove) eraseIndices(_state.variableNames, _columnsToRemove); } -// TODO move this into a simplifier class - -/// Turn constraints of the form ax <= b into an upper bound on x. -/// @returns false if the system is infeasible. -bool extractDirectConstraints(SolvingState& _state, bool& _changed) -{ - vector constraintsToRemove(_state.constraints.size(), false); - bool needsRemoval = false; - for (auto const& [index, constraint]: _state.constraints | ranges::views::enumerate) - { - auto nonzeroCoefficients = constraint.data | ranges::views::enumerate | ranges::views::tail | ranges::views::filter( - [](std::pair const& _x) { return !!_x.second; } - ); - // TODO we can exit early on in the loop above since we only care about zero, one or more than one nonzero entries. - // TODO could also use iterators and exit if we can advance it twice. - auto numNonzero = ranges::distance(nonzeroCoefficients); - if (numNonzero > 1) - continue; - constraintsToRemove[index] = true; - needsRemoval = true; - if (numNonzero == 0) - { - // 0 <= b or 0 = b - if ( - constraint.data.front().numerator() < 0 || - (constraint.equality && constraint.data.front()) - ) - return false; // Infeasible. - } - else - { - auto&& [varIndex, factor] = nonzeroCoefficients.front(); - // a * x <= b - rational bound = constraint.data[0] / factor; - if ( - (factor >= 0 || constraint.equality) && - (!_state.bounds[varIndex].upper || bound < _state.bounds[varIndex].upper) - ) - _state.bounds[varIndex].upper = bound; - if ( - (factor <= 0 || constraint.equality) && - (!_state.bounds[varIndex].lower || bound > _state.bounds[varIndex].lower) - ) - // Lower bound must be at least zero. - _state.bounds[varIndex].lower = max(rational{}, bound); - } - } - if (needsRemoval) - { - _changed = true; - eraseIndices(_state.constraints, constraintsToRemove); - } - return true; -} - -/// Remove variables that have equal lower and upper bound. -/// @returns false if the system is infeasible. -bool removeFixedVariables(SolvingState& _state, map& _model, bool& _changed) -{ - for (auto const& [index, bounds]: _state.bounds | ranges::views::enumerate) - { - if (!bounds.upper || (!bounds.lower && bounds.upper->numerator() > 0)) - continue; - // Lower bound must be at least zero. - rational lower = max(rational{}, bounds.lower ? *bounds.lower : rational{}); - rational upper = *bounds.upper; - if (upper < lower) - return false; // Infeasible. - if (upper != lower) - continue; - _model[_state.variableNames.at(index)] = lower; - _state.bounds[index] = {}; - _changed = true; - - // substitute variable - for (Constraint& constraint: _state.constraints) - if (constraint.data[index]) - { - constraint.data[0] -= constraint.data[index] * lower; - constraint.data[index] = 0; - } - } - - return true; -} - -bool removeEmptyColumns(SolvingState& _state, map& _model, bool& _changed) -{ - vector variablesSeen(_state.bounds.size(), false); - for (auto const& constraint: _state.constraints) - { - for (auto&& [index, factor]: constraint.data | ranges::views::enumerate | ranges::views::tail) - if (factor) - variablesSeen[index] = true; - } - - // TODO we could assert that any variable we remove does not have conflicting bounds. - // (We also remove the bounds). - - vector variablesToRemove(variablesSeen.size(), false); - bool needsRemoval = false; - for (auto&& [i, seen]: variablesSeen | ranges::views::enumerate | ranges::views::tail) - if (!seen) - { - variablesToRemove[i] = true; - needsRemoval = true; - // TODO actually it is unbounded if _state.bounds.at(i).upper is nullopt. - if (_state.bounds.at(i).lower || _state.bounds.at(i).upper) - _model[_state.variableNames.at(i)] = - _state.bounds.at(i).upper ? - *_state.bounds.at(i).upper : - *_state.bounds.at(i).lower; - } - if (needsRemoval) - { - _changed = true; - removeColumns(_state, variablesToRemove); - } - return true; -} - auto nonZeroEntriesInColumn(SolvingState const& _state, size_t _column) { return _state.constraints | ranges::views::enumerate | - ranges::views::filter([=](auto const& _entry) { return _entry.second.data[_column] != 0; }) | + ranges::views::filter([=](auto const& _entry) { return _entry.second.data[_column]; }) | ranges::views::transform([](auto const& _entry) { return _entry.first; }); } @@ -615,41 +494,6 @@ struct ProblemSplitter vector seenColumns; }; - -/// Simplifies the solving state according to some rules (remove rows without variables, etc). -/// @returns false if the state is determined to be infeasible during this process. -bool simplifySolvingState(SolvingState& _state, map& _model) -{ - // - Constraints with exactly one nonzero coefficient represent "a x <= b" - // and thus are turned into bounds. - // - Constraints with zero nonzero coefficients are constant relations. - // If such a relation is false, answer "infeasible", otherwise remove the constraint. - // - Empty columns can be removed. - // - Variables with matching bounds can be removed from the problem by substitution. - - bool changed = true; - while (changed) - { - changed = false; - - if (!removeFixedVariables(_state, _model, changed)) - return false; - - if (!extractDirectConstraints(_state, changed)) - return false; - - if (!removeFixedVariables(_state, _model, changed)) - return false; - - if (!removeEmptyColumns(_state, _model, changed)) - return false; - } - - // TODO return the values selected for named variables in order to - // be used when returning the model. - return true; -} - void normalizeRowLengths(SolvingState& _state) { size_t vars = max(_state.variableNames.size(), _state.bounds.size()); @@ -750,15 +594,160 @@ string SolvingState::toString() const return result; } +pair> SolvingStateSimplifier::simplify() +{ + + do + { + m_changed = false; + if ( + !removeFixedVariables() || + !extractDirectConstraints() || + // Used twice on purpose + !removeFixedVariables() || + !removeEmptyColumns() + ) + return {LPResult::Infeasible, {}}; + } + while (m_changed); + + return {LPResult::Unknown, move(m_model)}; +} + +bool SolvingStateSimplifier::removeFixedVariables() +{ + for (auto const& [index, bounds]: m_state.bounds | ranges::views::enumerate) + { + if (!bounds.upper || (!bounds.lower && bounds.upper->numerator() > 0)) + continue; + // Lower bound must be at least zero. + rational lower = max(rational{}, bounds.lower ? *bounds.lower : rational{}); + rational upper = *bounds.upper; + if (upper < lower) + return false; // Infeasible. + if (upper != lower) + continue; + m_model[m_state.variableNames.at(index)] = lower; + m_state.bounds[index] = {}; + m_changed = true; + + // substitute variable + for (Constraint& constraint: m_state.constraints) + if (constraint.data[index]) + { + constraint.data[0] -= constraint.data[index] * lower; + constraint.data[index] = 0; + } + } + + return true; +} + +bool SolvingStateSimplifier::extractDirectConstraints() +{ + vector constraintsToRemove(m_state.constraints.size(), false); + bool needsRemoval = false; + for (auto const& [index, constraint]: m_state.constraints | ranges::views::enumerate) + { + auto nonzeroCoefficients = constraint.data | ranges::views::enumerate | ranges::views::tail | ranges::views::filter( + [](std::pair const& _x) { return !!_x.second; } + ); + // TODO we can exit early on in the loop above since we only care about zero, one or more than one nonzero entries. + // TODO could also use iterators and exit if we can advance it twice. + auto numNonzero = ranges::distance(nonzeroCoefficients); + if (numNonzero > 1) + continue; + constraintsToRemove[index] = true; + needsRemoval = true; + if (numNonzero == 0) + { + // 0 <= b or 0 = b + if ( + constraint.data.front().numerator() < 0 || + (constraint.equality && constraint.data.front()) + ) + return false; // Infeasible. + } + else + { + auto&& [varIndex, factor] = nonzeroCoefficients.front(); + // a * x <= b + rational bound = constraint.data[0] / factor; + if ( + (factor >= 0 || constraint.equality) && + (!m_state.bounds[varIndex].upper || bound < m_state.bounds[varIndex].upper) + ) + m_state.bounds[varIndex].upper = bound; + if ( + (factor <= 0 || constraint.equality) && + (!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); + } + } + if (needsRemoval) + { + m_changed = true; + eraseIndices(m_state.constraints, constraintsToRemove); + } + return true; +} + +bool SolvingStateSimplifier::removeEmptyColumns() +{ + vector variablesSeen(m_state.bounds.size(), false); + for (auto const& constraint: m_state.constraints) + { + for (auto&& [index, factor]: constraint.data | ranges::views::enumerate | ranges::views::tail) + if (factor) + variablesSeen[index] = true; + } + + vector variablesToRemove(variablesSeen.size(), false); + bool needsRemoval = false; + for (auto&& [i, seen]: variablesSeen | ranges::views::enumerate | ranges::views::tail) + if (!seen) + { + variablesToRemove[i] = true; + needsRemoval = true; + SolvingState::Bounds const& bounds = m_state.bounds.at(i); + // TODO actually it is unbounded if m_state.bounds.at(i).upper is nullopt. + if (bounds.lower || bounds.upper) + { + solAssert(!bounds.upper || bounds.upper >= 0); + if (bounds.lower && bounds.upper) + solAssert(*bounds.lower <= *bounds.upper); + m_model[m_state.variableNames.at(i)] = + bounds.upper ? + *bounds.upper : + *bounds.lower; + } + } + if (needsRemoval) + { + m_changed = true; + removeColumns(m_state, variablesToRemove); + } + return true; +} + pair> LPSolver::check(SolvingState _state) { normalizeRowLengths(_state); - map model; - - if (!simplifySolvingState(_state, model)) + auto&& [simplificationResult, model] = SolvingStateSimplifier{_state}.simplify(); + switch (simplificationResult) + { + case LPResult::Infeasible: return {LPResult::Infeasible, {}}; + case LPResult::Feasible: + case LPResult::Unbounded: + solAssert(false); + case LPResult::Unknown: + break; + } bool canOnlyBeUnknown = false; ProblemSplitter splitter(_state); @@ -811,3 +800,4 @@ pair> LPSolver::check(SolvingState _state) return {LPResult::Feasible, move(model)}; } + diff --git a/libsolutil/LP.h b/libsolutil/LP.h index e9abb2392..31a9862a2 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -70,9 +70,53 @@ struct SolvingState enum class LPResult { Unknown, - Unbounded, - Feasible, - Infeasible + Unbounded, ///< System has a solution, but it can have an arbitrary objective value. + Feasible, ///< System has a solution (it might be unbounded, though). + Infeasible ///< System does not have any solution. +}; + + +/** + * Applies several strategies to simplify a given solving state. + * During these simplifications, it can sometimes already be determined if the + * state is feasible or not. + * Since some variables can be fixed to specific values, it returns a + * (partial) model. + * + * - Constraints with exactly one nonzero coefficient represent "a x <= b" + * and thus are turned into bounds. + * - Constraints with zero nonzero coefficients are constant relations. + * If such a relation is false, answer "infeasible", otherwise remove the constraint. + * - Empty columns can be removed. + * - Variables with matching bounds can be removed from the problem by substitution. + * + * Holds a reference to the solving state that is modified during operation. + */ +class SolvingStateSimplifier +{ +public: + SolvingStateSimplifier(SolvingState& _state): + m_state(_state) {} + + std::pair> simplify(); + +private: + /// Remove variables that have equal lower and upper bound. + /// @returns false if the system is infeasible. + bool 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(); + + /// Removes all-zeros columns. + bool removeEmptyColumns(); + + /// Set to true by the strategies if they performed some changes. + bool m_changed = false; + + SolvingState& m_state; + std::map m_model; }; /**