diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index 292912726..7bff849a2 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -270,15 +270,16 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom) //cerr << "Combining\n" << m_subProblems.at(_combineFrom)->toString(); //cerr << "\ninto\n" << m_subProblems.at(_combineInto)->toString(); SubProblem& combineInto = unseal(_combineInto); - SubProblem const& combineFrom = *m_subProblems[_combineFrom]; + SubProblem& combineFrom = *m_subProblems[_combineFrom]; size_t varShift = combineInto.variables.size(); #ifdef SPARSE size_t rowShift = combineInto.factors.rows(); + size_t colShift = combineInto.factors.columns(); for (size_t row = 0; row < combineFrom.factors.rows(); row++) - for (auto&& [col, v]: combineFrom.factors.iterateRow(row)) - combineinto.factors.insert(row + rowShift, col + colShift, move(v)); + for (auto&& entry: combineFrom.factors.iterateRow(row)) + combineInto.factors.insert(entry.row + rowShift, entry.col + colShift, move(entry.value)); #else size_t rowShift = combineInto.factors.size(); size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size(); @@ -362,7 +363,11 @@ void LPSolver::addConstraintToSubProblem( size_t slackIndex = addNewVariableToSubProblem(_subProblem); // Name is only needed for printing //problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++); +#ifdef SPARSE + problem.basicVariables[slackIndex] = problem.factors.rows(); +#else problem.basicVariables[slackIndex] = problem.factors.size(); +#endif if (_constraint.kind == Constraint::EQUAL) { problem.variables[slackIndex].bounds.lower = _constraint.data[0]; @@ -375,6 +380,41 @@ void LPSolver::addConstraintToSubProblem( // TODO it is a basic var, so we don't add it, unless we use this for basic vars. //problem.variablesPotentiallyOutOfBounds.insert(slackIndex); +#ifdef SPARSE + // Compress the constraint, i.e. turn outer variable indices into + // inner variable indices. + RationalWithDelta valueForSlack; + size_t row = problem.factors.rows(); + // First, handle the basic variables. + LinearExpression basicVarNullifier; + for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail()) + if (entry) + { + size_t innerIndex = problem.varMapping.at(outerIndex); + if (problem.basicVariables.count(innerIndex)) + { + problem.factors.addMultipleOfRow( + problem.basicVariables[innerIndex], + row, + entry + ); + // TODO could remove + problem.factors.entry(row, innerIndex).value = 0; + } + } + for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail()) + if (entry) + { + size_t innerIndex = problem.varMapping.at(outerIndex); + if (!problem.basicVariables.count(innerIndex)) + problem.factors.entry(row, innerIndex).value += entry; + valueForSlack += problem.variables[innerIndex].value * entry; + } + + problem.factors.entry(row, slackIndex).value = -1; + + problem.basicVariables[slackIndex] = row; +#else // Compress the constraint, i.e. turn outer variable indices into // inner variable indices. RationalWithDelta valueForSlack; @@ -399,7 +439,9 @@ void LPSolver::addConstraintToSubProblem( compressedConstraint[slackIndex] = -1; problem.factors.emplace_back(move(compressedConstraint)); + problem.basicVariables[slackIndex] = problem.factors.size() - 1; +#endif problem.variables[slackIndex].value = valueForSlack; } @@ -414,8 +456,10 @@ size_t LPSolver::addNewVariableToSubProblem(size_t _subProblem) { SubProblem& problem = unseal(_subProblem); size_t index = problem.variables.size(); +#ifndef SPARSE for (LinearExpression& c: problem.factors) c.resize(index + 1); +#endif problem.variables.emplace_back(); return index; } @@ -519,12 +563,23 @@ string LPSolver::SubProblem::toString() const resultString += " "; resultString += " := " + v.value.toString() + "\n"; } +#ifdef SPARSE + for (size_t rowIndex = 0; rowIndex < factors.rows(); rowIndex++) +#else for (auto&& [rowIndex, row]: factors | ranges::views::enumerate) +#endif { string basicVarPrefix; string rowString; +#ifdef SPARSE + for (auto&& entry: const_cast(factors).iterateRow(rowIndex)) + { + rational const& f = entry.value; + size_t i = entry.col; +#else for (auto&& [i, f]: row.enumerate()) { +#endif if (basicVariables.count(i) && basicVariables.at(i) == rowIndex) { solAssert(f == -1); @@ -590,13 +645,24 @@ void LPSolver::SubProblem::update(size_t _varIndex, RationalWithDelta const& _va { RationalWithDelta delta = _value - variables[_varIndex].value; variables[_varIndex].value = _value; +#ifdef SPARSE + // TODO can we store that? + map basicVarForRow = invertMap(basicVariables); + for (auto&& entry: factors.iterateColumn(_varIndex)) + if (entry.value && basicVarForRow.count(entry.row)) + { + size_t j = basicVarForRow[entry.row]; + variables[j].value += delta * entry.value; + //variablesPotentiallyOutOfBounds.insert(j); + } +#else for (auto&& [j, row]: basicVariables) if (factors[row][_varIndex]) { variables[j].value += delta * factors[row][_varIndex]; //variablesPotentiallyOutOfBounds.insert(j); } - +#endif } optional LPSolver::SubProblem::firstConflictingBasicVariable() const @@ -619,9 +685,16 @@ optional LPSolver::SubProblem::firstReplacementVar( bool _increasing ) const { +#ifdef SPARSE + for (auto&& entry: const_cast(factors).iterateRow(basicVariables.at(_basicVarToReplace))) + { + size_t i = entry.col; + rational const& factor = entry.value; +#else LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)]; for (auto const& [i, factor]: basicVarEquation.enumerate()) { +#endif if (i == _basicVarToReplace || !factor) continue; bool positive = factor > 0; @@ -647,9 +720,16 @@ set LPSolver::SubProblem::reasonsForUnsat( else if (!_increasing && variables[_basicVarToReplace].bounds.upperReason) r.insert(*variables[_basicVarToReplace].bounds.upperReason); +#ifdef SPARSE + for (auto&& entry: const_cast(factors).iterateRow(basicVariables.at(_basicVarToReplace))) + { + size_t i = entry.col; + rational const& factor = entry.value; +#else LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)]; for (auto const& [i, factor]: basicVarEquation.enumerate()) { +#endif if (i == _basicVarToReplace || !factor) continue; bool positive = factor > 0; @@ -669,6 +749,16 @@ void LPSolver::SubProblem::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 = basicVariables[_old]; +#ifdef SPARSE + rational pivot = factors.entry(pivotRow, _new).value; + solAssert(pivot != 0, ""); + if (pivot != -1) + factors.multiplyRowByFactor(pivotRow, rational{-1} / pivot); + + for (auto&& entry: factors.iterateColumn(_new)) + if (entry.row != pivotRow) + factors.addMultipleOfRow(pivotRow, entry.row, entry.value); +#else LinearExpression& pivotRowData = factors[pivotRow]; rational pivot = pivotRowData[_new]; @@ -691,6 +781,7 @@ void LPSolver::SubProblem::pivot(size_t _old, size_t _new) for (size_t i = 0; i < factors.size(); ++i) if (i != pivotRow) subtractMultipleOfPivotRow(factors[i]); +#endif basicVariables.erase(_old); basicVariables[_new] = pivotRow; @@ -702,14 +793,30 @@ void LPSolver::SubProblem::pivotAndUpdate( size_t _newBasicVar ) { +#ifdef SPARSE + RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(_oldBasicVar, _newBasicVar).value; +#else RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors[basicVariables[_oldBasicVar]][_newBasicVar]; +#endif variables[_oldBasicVar].value = _newValue; variables[_newBasicVar].value += theta; +#ifdef SPARSE + // TODO can we store that? + map basicVarForRow = invertMap(basicVariables); + for (auto&& entry: factors.iterateColumn(_newBasicVar)) + if (basicVarForRow.count(entry.row)) + { + size_t i = basicVarForRow[entry.row]; + if (i != _oldBasicVar) + variables[i].value += theta * entry.value; + } +#else for (auto const& [i, row]: basicVariables) if (i != _oldBasicVar && factors[row][_newBasicVar]) variables[i].value += theta * factors[row][_newBasicVar]; +#endif pivot(_oldBasicVar, _newBasicVar); } diff --git a/libsolutil/LinearExpression.cpp b/libsolutil/LinearExpression.cpp index b8c61fcc5..dd75b9ba1 100644 --- a/libsolutil/LinearExpression.cpp +++ b/libsolutil/LinearExpression.cpp @@ -23,10 +23,16 @@ using namespace std; SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::begin() { - return SparseMatrixIterator( - m_isRow ? m_matrix.m_row_start[m_RowOrColumn] : m_matrix.m_col_start[m_RowOrColumn], - m_isRow - ); + if ( + (m_isRow && m_rowOrColumn >= m_matrix.m_row_start.size()) || + (!m_isRow && m_rowOrColumn >= m_matrix.m_col_start.size()) + ) + return SparseMatrixIterator(nullptr, m_isRow); + else + return SparseMatrixIterator( + m_isRow ? m_matrix.m_row_start[m_rowOrColumn] : m_matrix.m_col_start[m_rowOrColumn], + m_isRow + ); } SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::end() @@ -34,7 +40,7 @@ SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::end() return SparseMatrixIterator(nullptr, m_isRow); } -SparseMatrix::IteratorCombiner SparseMatrix::enumerateColumn(size_t _column) +SparseMatrix::IteratorCombiner SparseMatrix::iterateColumn(size_t _column) { return IteratorCombiner{ _column, @@ -43,7 +49,7 @@ SparseMatrix::IteratorCombiner SparseMatrix::enumerateColumn(size_t _column) }; } -SparseMatrix::IteratorCombiner SparseMatrix::enumerateRow(size_t _row) +SparseMatrix::IteratorCombiner SparseMatrix::iterateRow(size_t _row) { return IteratorCombiner{ _row, @@ -67,6 +73,11 @@ void SparseMatrix::multiplyRowByFactor(size_t _row, rational const& _factor) void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor) { + if (_targetRow >= m_row_start.size()) + { + m_row_start.resize(_targetRow + 1); + m_row_end.resize(_targetRow + 1); + } Entry* source = m_row_start[_sourceRow]; Entry* target = m_row_start[_targetRow]; @@ -90,6 +101,30 @@ void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, ration } } +SparseMatrix::Entry& SparseMatrix::entry(size_t _row, size_t _column) +{ + Entry* successor = entryOrSuccessorInRow(_row, _column); + if (successor && successor->col == _column) + return *successor; + else + return *prependInRow(successor, _row, _column, {}); +} + +void SparseMatrix::insert(size_t _row, size_t _column, rational _value) +{ + if (_column >= m_col_start.size()) + { + m_col_start.resize(_column + 1); + m_col_end.resize(_column + 1); + } + if (_row >= m_row_start.size()) + { + m_row_start.resize(_row + 1); + m_row_end.resize(_row + 1); + } + + prependInRow(entryOrSuccessorInRow(_row, _column), _row, _column, move(_value)); +} void SparseMatrix::appendRow(LinearExpression const& _entries) { @@ -103,6 +138,19 @@ void SparseMatrix::appendRow(LinearExpression const& _entries) } } +SparseMatrix::Entry* SparseMatrix::entryOrSuccessorInRow(size_t _row, size_t _column) +{ + Entry* successor = nullptr; + if (m_row_end[_row] && m_row_end[_row]->col >= _column) + { + successor = m_row_start[_row]; + // TODO could choose to search from end + while (successor && successor->col < _column) + successor = successor->next_in_row; + } + return successor; +} + void SparseMatrix::remove(SparseMatrix::Entry& _e) { if (_e.prev_in_row) diff --git a/libsolutil/LinearExpression.h b/libsolutil/LinearExpression.h index 9c57f0c25..973804074 100644 --- a/libsolutil/LinearExpression.h +++ b/libsolutil/LinearExpression.h @@ -213,7 +213,7 @@ public: pointer operator->() { return m_ptr; } SparseMatrixIterator& operator++() { - m_ptr = m_isRow ? m_ptr->next_in_row : m_pt->next_in_col; + m_ptr = m_isRow ? m_ptr->next_in_row : m_ptr->next_in_col; return *this; } SparseMatrixIterator operator++(int) { SparseMatrixIterator tmp = *this; ++(*this); return tmp; } @@ -232,7 +232,7 @@ public: }; struct IteratorCombiner { - size_t m_RowOrColumn; + size_t m_rowOrColumn; bool m_isRow; SparseMatrix& m_matrix; SparseMatrixIterator begin(); @@ -242,25 +242,30 @@ public: size_t rows() const { return m_row_start.size(); } size_t columns() const { return m_col_start.size(); } - /// @returns (i, v) for all non-zero v in the column _column - void enumerateColumn(size_t _column); - /// @returns (i, v) for all non-zero v in the row _row - void enumerateRow(size_t _row); + /// @returns Entry for all non-zero v in the column _column + IteratorCombiner iterateColumn(size_t _column); + /// @returns Entry for all non-zero v in the row _row + IteratorCombiner iterateRow(size_t _row); void multiplyRowByFactor(size_t _row, rational const& _factor); void addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor); - rational entry(size_t _row, size_t _column) const; + Entry& entry(size_t _row, size_t _column); + /// Inserts the value at the row/rolumn. + /// Assumes the entry does not exist yet. void insert(size_t _row, size_t _column, rational _value); void appendRow(LinearExpression const& _entries); private: + /// @returns the entry at the row/column if it exists or its successor in the row. + Entry* entryOrSuccessorInRow(size_t _row, size_t _column); void remove(Entry& _entry); /// Prepends a new entry before the given element or at end of row if nullptr. Entry* prependInRow(Entry* _successor, size_t _row, size_t _column, rational _value); void adjustColumnProperties(Entry& _entry); - std::vector> m_elements; + // TODO unique_ptr? + std::vector> m_elements; std::vector m_row_start; std::vector m_col_start; std::vector m_row_end;