diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index f964b029b..99a8c5365 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -50,6 +50,16 @@ using rational = boost::rational; namespace { +/// Disjunctively combined two vectors of bools. +inline std::vector& operator|=(std::vector& _x, std::vector const& _y) +{ + solAssert(_x.size() == _y.size(), ""); + for (size_t i = 0; i < _x.size(); ++i) + if (_y[i]) + _x[i] = true; + return _x; +} + /** * Simplex tableau. */ @@ -57,36 +67,39 @@ struct Tableau { /// The factors of the objective function (first row of the tableau) LinearExpression objective; - /// The tableau matrix (equational forrm). + /// The tableau matrix (equational form). std::vector data; }; /// Adds slack variables to remove non-equality costraints from a set of constraints. -/// The second return variable is true if new variables have been added. +/// The second return variable is true if a non-equality constraint was found and thus new variables have been added. pair, bool> toEquationalForm(vector _constraints) { size_t varsNeeded = static_cast(ranges::count_if(_constraints, [](Constraint const& _c) { return !_c.equality; })); if (varsNeeded > 0) { size_t columns = _constraints.at(0).data.size(); - size_t currentVariable = 0; + size_t varsAdded = 0; for (Constraint& constraint: _constraints) { solAssert(constraint.data.size() == columns, ""); - constraint.data.factors += vector(varsNeeded, rational{}); + constraint.data.resize(columns + varsNeeded); if (!constraint.equality) { constraint.equality = true; - constraint.data[columns + currentVariable] = bigint(1); - currentVariable++; + constraint.data[columns + varsAdded] = bigint(1); + varsAdded++; } } + solAssert(varsAdded == varsNeeded); } return make_pair(move(_constraints), varsNeeded > 0); } +/// Finds the simplex pivot column: The column with the largest positive objective factor. +/// If all objective factors are zero or negative, the optimum has been found and nullopt is returned. optional findPivotColumn(Tableau const& _tableau) { auto&& [maxColumn, maxValue] = ranges::max( @@ -101,6 +114,10 @@ optional findPivotColumn(Tableau const& _tableau) return maxColumn; } +/// Finds the simplex pivot row, given the column: +/// If there is no positive factor in the column, the problem is unbounded, nullopt is returned. +/// Otherwise, returns the row i such that c[i] / x[i] is minimal and x[i] is positive, where +/// c[i] is the constant factor (not the objective factor!) in row i. optional findPivotRow(Tableau const& _tableau, size_t _pivotColumn) { auto positiveColumnEntries = @@ -109,7 +126,7 @@ optional findPivotRow(Tableau const& _tableau, size_t _pivotColumn) return make_pair(i, _tableau.data[i][_pivotColumn]); }) | ranges::views::filter([](pair const& _entry) { - return _entry.second > 0; + return _entry.second.numerator() > 0; }); if (positiveColumnEntries.empty()) return nullopt; // unbounded @@ -124,7 +141,7 @@ optional findPivotRow(Tableau const& _tableau, size_t _pivotColumn) } /// Performs equivalence transform on @a _tableau, so that -/// the column @a _pivotColumn is all zeros except for @a _pivotRow, +/// the column @a _pivotColumn is all zeros (including the objective row) except for @a _pivotRow, /// where it is 1. void performPivot(Tableau& _tableau, size_t _pivotRow, size_t _pivotColumn) { @@ -135,23 +152,27 @@ void performPivot(Tableau& _tableau, size_t _pivotRow, size_t _pivotColumn) solAssert(_tableau.data[_pivotRow][_pivotColumn] == rational(1), ""); LinearExpression const& _pivotRowData = _tableau.data[_pivotRow]; - auto subtractPivotRow = [&](LinearExpression& _row) { + auto subtractMultipleOfPivotRow = [&](LinearExpression& _row) { if (_row[_pivotColumn] == rational{1}) _row -= _pivotRowData; - else if (_row[_pivotColumn] != rational{}) + else if (_row[_pivotColumn]) _row -= _row[_pivotColumn] * _pivotRowData; }; - subtractPivotRow(_tableau.objective); + subtractMultipleOfPivotRow(_tableau.objective); for (size_t i = 0; i < _tableau.data.size(); ++i) if (i != _pivotRow) - subtractPivotRow(_tableau.data[i]); + subtractMultipleOfPivotRow(_tableau.data[i]); } +/// Transforms the tableau such that the last vectors are basis vectors +/// and their objective coefficients are zero. +/// Makes various assumptions and should only be used after adding +/// a certain number of slack variables. void selectLastVectorsAsBasis(Tableau& _tableau) { // We might skip the operation for a column if it is already the correct - // unit vector and its cost coefficient is zero. + // unit vector and its objective coefficient is zero. size_t columns = _tableau.objective.size(); size_t rows = _tableau.data.size(); for (size_t i = 0; i < rows; ++i) @@ -161,7 +182,7 @@ void selectLastVectorsAsBasis(Tableau& _tableau) /// If column @a _column inside tableau is a basis vector /// (i.e. one entry is 1, the others are 0), returns the index /// of the 1, otherwise nullopt. -optional basisVariable(Tableau const& _tableau, size_t _column) +optional basisIndex(Tableau const& _tableau, size_t _column) { optional row; for (size_t i = 0; i < _tableau.data.size(); ++i) @@ -186,7 +207,7 @@ vector solutionVector(Tableau const& _tableau) vector rowsSeen(_tableau.data.size(), false); for (size_t j = 1; j < _tableau.objective.size(); j++) { - optional row = basisVariable(_tableau, j); + optional row = basisIndex(_tableau, j); if (row && rowsSeen[*row]) row = nullopt; result.emplace_back(row ? _tableau.data[*row][0] : rational{}); @@ -201,6 +222,7 @@ vector solutionVector(Tableau const& _tableau) /// Here, c is _tableau.objective and the first column of _tableau.data /// encodes b and the other columns encode A /// Assumes the tableau has a trivial basic feasible solution. +/// Tries for a number of iterations and then gives up. pair simplexEq(Tableau _tableau) { size_t const iterations = min(60, 50 + _tableau.objective.size() * 2); @@ -237,7 +259,7 @@ pair simplexPhaseI(Tableau _tableau) { if (_tableau.data[i][0] < 0) _tableau.data[i] *= -1; - _tableau.data[i].factors += vector(rows, bigint{}); + _tableau.data[i].resize(columns + rows); _tableau.data[i][columns + i] = 1; } _tableau.objective.factors = diff --git a/libsolutil/LP.h b/libsolutil/LP.h index 70df74736..e9abb2392 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -48,8 +48,8 @@ struct Constraint */ struct SolvingState { - /// Names of variables, the index zero should be left empty - /// (because zero corresponds to constants). + /// Names of variables. The index zero should be left empty + /// because zero corresponds to constants. std::vector variableNames; struct Bounds { diff --git a/libsolutil/LinearExpression.h b/libsolutil/LinearExpression.h index d2ed920ed..5002c82aa 100644 --- a/libsolutil/LinearExpression.h +++ b/libsolutil/LinearExpression.h @@ -39,7 +39,7 @@ using rational = boost::rational; /** * A linear expression of the form * factors[0] + factors[1] * X1 + factors[2] * X2 + ... - * where the variables X_i are implicit. + * The set and order of variables is implied. */ struct LinearExpression { @@ -82,6 +82,7 @@ struct LinearExpression factors.resize(_size); } + /// Sets the factor at @a _index to @a _factor and enlarges if needed. void resizeAndSet(size_t _index, rational _factor) { if (factors.size() <= _index) @@ -89,6 +90,7 @@ struct LinearExpression factors[_index] = move(_factor); } + /// @returns true if all factors of variables are zero. bool isConstant() const { return ranges::all_of(factors | ranges::views::tail, [](rational const& _v) { return !_v; }); @@ -153,10 +155,9 @@ struct LinearExpression } - /// Multiply two vectors where the first element of each vector is a constant factor. - /// Only works if at most one of the vector has a nonzero element after the first. - /// If this condition is violated, returns nullopt. - static std::optional vectorProduct( + /// Multiply two linear expression. This only works if at least one of them is a constant. + /// Returns nullopt otherwise. + friend std::optional operator*( std::optional _x, std::optional _y ) @@ -178,16 +179,4 @@ struct LinearExpression std::vector factors; }; -// TODO - - -inline std::vector& operator|=(std::vector& _x, std::vector const& _y) -{ - solAssert(_x.size() == _y.size(), ""); - for (size_t i = 0; i < _x.size(); ++i) - if (_y[i]) - _x[i] = true; - return _x; -} - }