Turn into class.

This commit is contained in:
chriseth 2022-04-01 19:01:50 +02:00
parent e2eeed6af1
commit 2ef255bdab
2 changed files with 219 additions and 221 deletions

View File

@ -381,225 +381,6 @@ pair<LPResult, vector<rational>> simplex(vector<Constraint> _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<size_t, size_t> normalizeForSolvingWithBounds(SolvingState& _state)
{
size_t varsNeeded = _state.constraints.size();
/*
for (Constraint& c: _state.constraints)
if (!c.equality || c.data[0])
varsNeeded++;
*/
map<size_t, size_t> 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<rational>& _assignments,
map<size_t, size_t> 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<size_t> firstConflictingBasicVariable(
SolvingState& _state,
vector<rational>& _assignments,
map<size_t, size_t> 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<size_t> firstReplacementVar(
SolvingState& _state,
vector<rational>& _assignments,
map<size_t, size_t> 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<size_t, size_t>& _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<size_t, size_t>& _basicVariables,
vector<rational>& _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<LPResult, vector<rational>> simplexWithBounds(SolvingState _state)
{
cout << "===========================" << endl;
cout << "===========================" << endl;
cout << "simpl with bounds on\n" << _state.toString() << endl;
map<size_t, size_t> basicVariables = normalizeForSolvingWithBounds(_state);
cout << "After norm ------------------------\n" << _state.toString() << endl;
vector<rational> 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<rational>{});
}
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<rational>{});
}
}
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<ReasonSet> 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<size_t> 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<size_t> 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<LPResult, variant<map<size_t, rational>, ReasonSet>> SolvingStateSimplifier::simplify()
{
do
@ -1150,8 +1124,7 @@ pair<LPResult, variant<Model, ReasonSet>> 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));
}

View File

@ -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<size_t> firstConflictingBasicVariable() const;
std::optional<size_t> 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<rational> m_assignments;
/// Variable index to row it controls.
std::map<size_t, size_t> m_basicVariables;
};
/**
* Applies several strategies to simplify a given solving state.
* During these simplifications, it can sometimes already be determined if the