From bba2832e5b45fc1cb7847a5b30faa3b33d09e9ae Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 11 Sep 2017 22:17:42 +0200 Subject: [PATCH] finished Walker-Chae method --- .../solver/NativeLinearEquationSolver.cpp | 113 +++++++++--------- src/storm/solver/NativeLinearEquationSolver.h | 20 +++- src/storm/storage/SparseMatrix.cpp | 42 +++++-- src/storm/storage/SparseMatrix.h | 11 +- .../NativeDtmcPrctlModelCheckerTest.cpp | 62 +++++----- 5 files changed, 139 insertions(+), 109 deletions(-) diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index b3c1d4f37..855b56fed 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -206,12 +206,19 @@ namespace storm { return converged; } - + + template + NativeLinearEquationSolver::WalkerChaeData::WalkerChaeData(storm::storage::SparseMatrix const& originalMatrix, std::vector const& originalB) : t(storm::utility::convertNumber(1000.0)) { + computeWalkerChaeMatrix(originalMatrix); + computeNewB(originalB); + precomputeAuxiliaryData(); + } + template - void NativeLinearEquationSolver::computeWalkerChaeMatrix() const { - storm::storage::BitVector columnsWithNegativeEntries(this->A->getColumnCount()); + void NativeLinearEquationSolver::WalkerChaeData::computeWalkerChaeMatrix(storm::storage::SparseMatrix const& originalMatrix) { + storm::storage::BitVector columnsWithNegativeEntries(originalMatrix.getColumnCount()); ValueType zero = storm::utility::zero(); - for (auto const& e : *this->A) { + for (auto const& e : originalMatrix) { if (e.getValue() < zero) { columnsWithNegativeEntries.set(e.getColumn()); } @@ -222,10 +229,10 @@ namespace storm { storm::storage::SparseMatrixBuilder builder; uint64_t row = 0; - for (; row < this->A->getRowCount(); ++row) { - for (auto const& entry : this->A->getRow(row)) { + for (; row < originalMatrix.getRowCount(); ++row) { + for (auto const& entry : originalMatrix.getRow(row)) { if (entry.getValue() < zero) { - builder.addNextValue(row, this->A->getRowCount() + columnsWithNegativeEntriesBefore[entry.getColumn()], -entry.getValue()); + builder.addNextValue(row, originalMatrix.getRowCount() + columnsWithNegativeEntriesBefore[entry.getColumn()], -entry.getValue()); } else { builder.addNextValue(row, entry.getColumn(), entry.getValue()); } @@ -234,11 +241,27 @@ namespace storm { ValueType one = storm::utility::one(); for (auto column : columnsWithNegativeEntries) { builder.addNextValue(row, column, one); - builder.addNextValue(row, this->A->getRowCount() + columnsWithNegativeEntriesBefore[column], one); + builder.addNextValue(row, originalMatrix.getRowCount() + columnsWithNegativeEntriesBefore[column], one); ++row; } - walkerChaeMatrix = std::make_unique>(builder.build()); + matrix = builder.build(); + } + + template + void NativeLinearEquationSolver::WalkerChaeData::computeNewB(std::vector const& originalB) { + b = std::vector(originalB); + b.resize(matrix.getRowCount()); + } + + template + void NativeLinearEquationSolver::WalkerChaeData::precomputeAuxiliaryData() { + columnSums = std::vector(matrix.getColumnCount()); + for (auto const& e : matrix) { + columnSums[e.getColumn()] += e.getValue(); + } + + newX.resize(matrix.getRowCount()); } template @@ -246,69 +269,45 @@ namespace storm { STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (WalkerChae)"); // (1) Compute an equivalent equation system that has only non-negative coefficients. - if (!walkerChaeMatrix) { - std::cout << *this->A << std::endl; - computeWalkerChaeMatrix(); - std::cout << *walkerChaeMatrix << std::endl; + if (!walkerChaeData) { + walkerChaeData = std::make_unique(*this->A, b); } // (2) Enlarge the vectors x and b to account for additional variables. - x.resize(walkerChaeMatrix->getRowCount()); - - if (walkerChaeMatrix->getRowCount() > this->A->getRowCount() && !walkerChaeB) { - walkerChaeB = std::make_unique>(b); - walkerChaeB->resize(x.size()); - } - - // Choose a value for t in the algorithm. - ValueType t = storm::utility::convertNumber(1000); - - // Precompute some data. - std::vector columnSums(x.size()); - for (auto const& e : *walkerChaeMatrix) { - STORM_LOG_ASSERT(e.getValue() >= storm::utility::zero(), "Expecting only non-negative entries in WalkerChae matrix."); - columnSums[e.getColumn()] += e.getValue(); - } + x.resize(walkerChaeData->matrix.getRowCount()); // Square the error bound, so we can use it to check for convergence. We take the squared error, because we // do not want to compute the root in the 2-norm computation. ValueType squaredErrorBound = storm::utility::pow(this->getSettings().getPrecision(), 2); - // Create a vector that always holds Ax. - std::vector currentAx(x.size()); - walkerChaeMatrix->multiplyWithVector(x, currentAx); - - // Create an auxiliary vector that intermediately stores the result of the Walker-Chae step. - std::vector tmpX(x.size()); - // Set up references to the x-vectors used in the iteration loop. std::vector* currentX = &x; - std::vector* nextX = &tmpX; + std::vector* nextX = &walkerChaeData->newX; + + std::vector tmp = walkerChaeData->matrix.getRowSumVector(); + storm::utility::vector::applyPointwise(tmp, walkerChaeData->b, walkerChaeData->b, [this] (ValueType const& first, ValueType const& second) { return walkerChaeData->t * first + second; } ); + + // Add t to all entries of x. + storm::utility::vector::applyPointwise(x, x, [this] (ValueType const& value) { return value + walkerChaeData->t; }); - // Prepare a function that adds t to its input. - auto addT = [t] (ValueType const& value) { return value + t; }; + // Create a vector that always holds Ax. + std::vector currentAx(x.size()); + walkerChaeData->matrix.multiplyWithVector(*currentX, currentAx); // (3) Perform iterations until convergence. bool converged = false; uint64_t iterations = 0; while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations()) { // Perform one Walker-Chae step. - A->performWalkerChaeStep(*currentX, columnSums, *walkerChaeB, currentAx, *nextX); - + walkerChaeData->matrix.performWalkerChaeStep(*currentX, walkerChaeData->columnSums, walkerChaeData->b, currentAx, *nextX); + // Compute new Ax. - A->multiplyWithVector(*nextX, currentAx); + walkerChaeData->matrix.multiplyWithVector(*nextX, currentAx); // Check for convergence. - converged = storm::utility::vector::computeSquaredNorm2Difference(currentAx, *walkerChaeB); - - // If the method did not yet converge, we need to update the value of Ax. - if (!converged) { - // TODO: scale matrix diagonal entries with t and add them to *walkerChaeB. - - // Add t to all entries of x. - storm::utility::vector::applyPointwise(x, x, addT); - } + converged = storm::utility::vector::computeSquaredNorm2Difference(currentAx, walkerChaeData->b) <= squaredErrorBound; + // Swap the x vectors for the next iteration. std::swap(currentX, nextX); // Increase iteration count so we can abort if convergence is too slow. @@ -317,14 +316,10 @@ namespace storm { // If the last iteration did not write to the original x we have to swap the contents, because the // output has to be written to the input parameter x. - if (currentX == &tmpX) { + if (currentX == &walkerChaeData->newX) { std::swap(x, *currentX); } - - if (!this->isCachingEnabled()) { - clearCache(); - } - + if (converged) { STORM_LOG_INFO("Iterative solver converged in " << iterations << " iterations."); } else { @@ -335,7 +330,7 @@ namespace storm { x.resize(this->A->getRowCount()); // Finalize solution vector. - storm::utility::vector::applyPointwise(x, x, [&t,iterations] (ValueType const& value) { return value - iterations * t; } ); + storm::utility::vector::applyPointwise(x, x, [this] (ValueType const& value) { return value - walkerChaeData->t; } ); if (!this->isCachingEnabled()) { clearCache(); @@ -424,7 +419,7 @@ namespace storm { template void NativeLinearEquationSolver::clearCache() const { jacobiDecomposition.reset(); - walkerChaeMatrix.reset(); + walkerChaeData.reset(); LinearEquationSolver::clearCache(); } diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index ebe180389..2c9633c55 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -68,8 +68,6 @@ namespace storm { virtual bool solveEquationsSOR(std::vector& x, std::vector const& b, ValueType const& omega) const; virtual bool solveEquationsJacobi(std::vector& x, std::vector const& b) const; virtual bool solveEquationsWalkerChae(std::vector& x, std::vector const& b) const; - - void computeWalkerChaeMatrix() const; // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. @@ -85,8 +83,22 @@ namespace storm { // cached auxiliary data mutable std::unique_ptr, std::vector>> jacobiDecomposition; - mutable std::unique_ptr> walkerChaeMatrix; - mutable std::unique_ptr> walkerChaeB; + struct WalkerChaeData { + WalkerChaeData(storm::storage::SparseMatrix const& originalMatrix, std::vector const& originalB); + + void computeWalkerChaeMatrix(storm::storage::SparseMatrix const& originalMatrix); + void computeNewB(std::vector const& originalB); + void precomputeAuxiliaryData(); + + storm::storage::SparseMatrix matrix; + std::vector b; + ValueType t; + + // Auxiliary data. + std::vector columnSums; + std::vector newX; + }; + mutable std::unique_ptr walkerChaeData; }; template diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 812e1214f..c622ed7d9 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -781,6 +781,18 @@ namespace storm { } } + template + std::vector SparseMatrix::getRowSumVector() const { + std::vector result(this->getRowCount()); + + index_type row = 0; + for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++row) { + *resultIt = getRowSum(row); + } + + return result; + } + template ValueType SparseMatrix::getConstrainedRowSum(index_type row, storm::storage::BitVector const& constraint) const { ValueType result = storm::utility::zero(); @@ -1465,25 +1477,31 @@ namespace storm { const_iterator it = this->begin(); const_iterator ite; std::vector::const_iterator rowIterator = rowIndications.begin(); - typename std::vector::iterator resultIterator = result.begin(); - typename std::vector::iterator resultIteratorEnd = result.end(); - typename std::vector::const_iterator xIterator = x.begin(); - typename std::vector::const_iterator columnSumsIterator = columnSums.begin(); - - for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++xIterator, ++columnSumsIterator) { - *resultIterator = storm::utility::zero(); - + + // Clear all previous entries. + ValueType zero = storm::utility::zero(); + for (auto& entry : result) { + entry = zero; + } + + for (index_type row = 0; row < rowCount; ++row, ++rowIterator) { for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { - *resultIterator += it->getValue() * (b[it->getColumn()] / ax[it->getColumn()]); + result[it->getColumn()] += it->getValue() * (b[row] / ax[row]); } - - *resultIterator *= (*xIterator / *columnSumsIterator); + } + + auto xIterator = x.begin(); + auto sumsIterator = columnSums.begin(); + for (auto& entry : result) { + entry *= *xIterator / *sumsIterator; + ++xIterator; + ++sumsIterator; } } #ifdef STORM_HAVE_CARL template<> - void SparseMatrix::performWalkerChaeStep(std::vector const& x, std::vector const& columnSums, std::vector const& b, std::vector const& ax, std::vector& result) const { + void SparseMatrix::performWalkerChaeStep(std::vector const& x, std::vector const& rowSums, std::vector const& b, std::vector const& ax, std::vector& result) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index e5eb9ae9f..bea7602d8 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -608,6 +608,13 @@ namespace storm { */ void makeRowDirac(index_type row, index_type column); + /* + * Sums the entries in all rows. + * + * @return The vector of sums of the entries in the respective rows. + */ + std::vector getRowSumVector() const; + /* * Sums the entries in the given row and columns. * @@ -838,7 +845,7 @@ namespace storm { * @param divisors The divisors with which each row is divided. */ void divideRowsInPlace(std::vector const& divisors); - + /*! * Performs one step of the successive over-relaxation technique. * @@ -852,7 +859,7 @@ namespace storm { * Performs one step of the Walker-Chae technique. * * @param x The current solution vector. - * @param columnSums The sums of the entries of the individual columns. + * @param columnSums The sums the individual columns. * @param b The 'right-hand side' of the problem. * @param ax A vector resulting from multiplying the current matrix with the vector x. * @param result The vector to which to write the result. diff --git a/src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp index ab27a9b01..24e72c778 100644 --- a/src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp +++ b/src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp @@ -194,38 +194,35 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) { EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule().getPrecision()); } -// Does not converge any more. :( -// { -// matrixBuilder = storm::storage::SparseMatrixBuilder(3, 3, 3); -// matrixBuilder.addNextValue(0, 1, 1); -// matrixBuilder.addNextValue(1, 2, 1); -// matrixBuilder.addNextValue(2, 0, 1); -// storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); -// -// storm::models::sparse::StateLabeling ap(3); -// ap.addLabel("a"); -// ap.addLabelToState("a", 2); -// -// dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); -// -// auto factory = std::make_unique>(); -// factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings::SolutionMethod::SOR); -// factory->getSettings().setOmega(0.9); -// storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::move(factory)); -// -// std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); -// -// std::unique_ptr result = checker.check(*formula); -// storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); -// -// EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule().getPrecision()); -// EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule().getPrecision()); -// EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule().getPrecision()); -// } + { + matrixBuilder = storm::storage::SparseMatrixBuilder(3, 3, 3); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 2, 1); + matrixBuilder.addNextValue(2, 0, 1); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 2); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); + + auto factory = std::make_unique>(); + factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings::SolutionMethod::WalkerChae); + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::move(factory)); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule().getPrecision()); + } } -// Test is disabled as it does not converge any more. :( -TEST(DISABLED_NativeDtmcPrctlModelCheckerTest, LRA) { +TEST(NativeDtmcPrctlModelCheckerTest, LRA) { storm::storage::SparseMatrixBuilder matrixBuilder; std::shared_ptr> dtmc; @@ -276,8 +273,9 @@ TEST(DISABLED_NativeDtmcPrctlModelCheckerTest, LRA) { dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); auto factory = std::make_unique>(); - factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings::SolutionMethod::SOR); - factory->getSettings().setOmega(0.9); + factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings::SolutionMethod::WalkerChae); + factory->getSettings().setPrecision(1e-7); + factory->getSettings().setMaximalNumberOfIterations(50000); storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::move(factory)); std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");