From 594130dd50671c2b348bf455de4483212b32f359 Mon Sep 17 00:00:00 2001 From: chriseth Date: Thu, 18 Aug 2022 19:14:42 +0200 Subject: [PATCH] fixes --- libsolutil/LP.cpp | 28 +++++++++---- libsolutil/LP.h | 1 + libsolutil/LinearExpression.cpp | 73 +++++++++++++++------------------ libsolutil/LinearExpression.h | 7 +--- 4 files changed, 58 insertions(+), 51 deletions(-) diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index 7bff849a2..a51592fcd 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -51,6 +51,7 @@ using namespace solidity::util; using rational = boost::rational; +#define DEBUG namespace { @@ -279,7 +280,7 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom) for (size_t row = 0; row < combineFrom.factors.rows(); row++) for (auto&& entry: combineFrom.factors.iterateRow(row)) - combineInto.factors.insert(entry.row + rowShift, entry.col + colShift, move(entry.value)); + combineInto.factors.entry(entry.row + rowShift, entry.col + colShift).value = move(entry.value); #else size_t rowShift = combineInto.factors.size(); size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size(); @@ -362,7 +363,9 @@ void LPSolver::addConstraintToSubProblem( // Introduce the slack variable. size_t slackIndex = addNewVariableToSubProblem(_subProblem); // Name is only needed for printing - //problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++); +#ifdef DEBUG + problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++); +#endif #ifdef SPARSE problem.basicVariables[slackIndex] = problem.factors.rows(); #else @@ -398,8 +401,7 @@ void LPSolver::addConstraintToSubProblem( row, entry ); - // TODO could remove - problem.factors.entry(row, innerIndex).value = 0; + problem.factors.remove(problem.factors.entry(row, innerIndex)); } } for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail()) @@ -407,7 +409,12 @@ void LPSolver::addConstraintToSubProblem( { size_t innerIndex = problem.varMapping.at(outerIndex); if (!problem.basicVariables.count(innerIndex)) - problem.factors.entry(row, innerIndex).value += entry; + { + SparseMatrix::Entry& e = problem.factors.entry(row, innerIndex); + e.value += entry; + if (!e.value) + problem.factors.remove(e); + } valueForSlack += problem.variables[innerIndex].value * entry; } @@ -506,6 +513,7 @@ LPResult LPSolver::SubProblem::check() { #ifdef DEBUG cerr << "Replacing by " << variables[*replacementVar].name << endl; + cerr << "Setting basic var to to " << basicVar.bounds.lower->m_main << endl; #endif pivotAndUpdate(*bvi, *basicVar.bounds.lower, *replacementVar); @@ -575,6 +583,7 @@ string LPSolver::SubProblem::toString() const for (auto&& entry: const_cast(factors).iterateRow(rowIndex)) { rational const& f = entry.value; + solAssert(!!f); size_t i = entry.col; #else for (auto&& [i, f]: row.enumerate()) @@ -755,9 +764,14 @@ void LPSolver::SubProblem::pivot(size_t _old, size_t _new) if (pivot != -1) factors.multiplyRowByFactor(pivotRow, rational{-1} / pivot); - for (auto&& entry: factors.iterateColumn(_new)) + for (auto it = factors.iterateColumn(_new).begin(); it != factors.iterateColumn(_new).end(); ) + { + SparseMatrix::Entry& entry = *it; + // Increment becasue "addMultipleOfRow" might invalidate the iterator + ++it; if (entry.row != pivotRow) factors.addMultipleOfRow(pivotRow, entry.row, entry.value); + } #else LinearExpression& pivotRowData = factors[pivotRow]; @@ -794,7 +808,7 @@ void LPSolver::SubProblem::pivotAndUpdate( ) { #ifdef SPARSE - RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(_oldBasicVar, _newBasicVar).value; + RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(basicVariables[_oldBasicVar], _newBasicVar).value; #else RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors[basicVariables[_oldBasicVar]][_newBasicVar]; #endif diff --git a/libsolutil/LP.h b/libsolutil/LP.h index df5ba7da2..bed034243 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -19,6 +19,7 @@ // use sparse matrices #define SPARSE 1 +#define DEBUG #include #include diff --git a/libsolutil/LinearExpression.cpp b/libsolutil/LinearExpression.cpp index 412b1d23e..e947b807b 100644 --- a/libsolutil/LinearExpression.cpp +++ b/libsolutil/LinearExpression.cpp @@ -71,13 +71,10 @@ void SparseMatrix::multiplyRowByFactor(size_t _row, rational const& _factor) } } -void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, 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); - } + ensureSize(_targetRow, 0); + solAssert(_sourceRow != _targetRow); Entry* source = m_row_start[_sourceRow]; Entry* target = m_row_start[_targetRow]; @@ -94,8 +91,13 @@ void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, ration remove(*target); target = next; } + else + target = target->next_in_row; } - target = prependInRow(target, _targetRow, source->col, _factor * source->value)->next_in_row; + else if (rational newValue = _factor * source->value) + target = prependInRow(target, _targetRow, source->col, newValue)->next_in_row; + else if (target) + target = target->next_in_row; source = source->next_in_row; } @@ -111,18 +113,35 @@ SparseMatrix::Entry& SparseMatrix::entry(size_t _row, size_t _column) return *prependInRow(successor, _row, _column, {}); } -void SparseMatrix::insert(size_t _row, size_t _column, rational _value) +void SparseMatrix::remove(SparseMatrix::Entry& _e) { - ensureSize(_row, _column); - prependInRow(entryOrSuccessorInRow(_row, _column), _row, _column, move(_value)); + // TODO this does not deallocate the entry. + // At some point we should perform cleanup + // and swap entries in the vector + if (_e.prev_in_row) + _e.prev_in_row->next_in_row = _e.next_in_row; + else + m_row_start[_e.row] = _e.next_in_row; + if (_e.next_in_row) + _e.next_in_row->prev_in_row = _e.prev_in_row; + else + m_row_end[_e.row] = _e.prev_in_row; + if (_e.prev_in_col) + _e.prev_in_col->next_in_col = _e.next_in_col; + else + m_col_start[_e.col] = _e.next_in_col; + if (_e.next_in_col) + _e.next_in_col->prev_in_col = _e.prev_in_col; + else + m_col_end[_e.col] = _e.prev_in_col; } void SparseMatrix::appendRow(LinearExpression const& _entries) { - m_row_start.push_back(nullptr); - m_row_end.push_back(nullptr); - size_t row_nr = m_row_start.size() - 1; - for (auto&& [i, v]: _entries.enumerate()) { + size_t row_nr = rows(); + ensureSize(row_nr, 0); + for (auto&& [i, v]: _entries.enumerate()) + { if (!v) continue; prependInRow(nullptr, row_nr, i, move(v)); @@ -156,26 +175,6 @@ SparseMatrix::Entry* SparseMatrix::entryOrSuccessorInRow(size_t _row, size_t _co return successor; } -void SparseMatrix::remove(SparseMatrix::Entry& _e) -{ - if (_e.prev_in_row) - _e.prev_in_row->next_in_row = _e.next_in_row; - else - m_row_start[_e.row] = _e.next_in_row; - if (_e.next_in_row) - _e.next_in_row->prev_in_row = _e.prev_in_row; - else - m_row_end[_e.row] = _e.prev_in_row; - if (_e.prev_in_col) - _e.prev_in_col->next_in_col = _e.next_in_col; - else - m_col_start[_e.col] = _e.next_in_col; - if (_e.next_in_col) - _e.next_in_col->prev_in_col = _e.prev_in_col; - else - m_col_end[_e.col] = _e.prev_in_col; -} - SparseMatrix::Entry* SparseMatrix::prependInRow(Entry* _successor, size_t _row, size_t _column, rational _value) { m_elements.emplace_back(make_unique(Entry{ @@ -211,11 +210,7 @@ void SparseMatrix::adjustColumnProperties(Entry& _entry) { size_t column = _entry.col; - if (column >= m_col_start.size()) - { - m_col_start.resize(column + 1); - m_col_end.resize(column + 1); - } + solAssert(column < m_col_start.size()); Entry* c = nullptr; if (m_col_end[column] && m_col_end[column]->row > _entry.row) { diff --git a/libsolutil/LinearExpression.h b/libsolutil/LinearExpression.h index 06b0bf157..cdf8f599a 100644 --- a/libsolutil/LinearExpression.h +++ b/libsolutil/LinearExpression.h @@ -247,11 +247,9 @@ public: /// @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); + void addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational _factor); 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 remove(Entry& _entry); void appendRow(LinearExpression const& _entries); @@ -260,7 +258,6 @@ 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);