sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

finished Walker-Chae method

main
dehnert 8 years ago
parent
5440d164b2
commit
bba2832e5b
5 changed files with 139 additions and 109 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 113
      src/storm/solver/NativeLinearEquationSolver.cpp
  2. 20
      src/storm/solver/NativeLinearEquationSolver.h
  3. 42
      src/storm/storage/SparseMatrix.cpp
  4. 11
      src/storm/storage/SparseMatrix.h
  5. 62
      src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp

113
src/storm/solver/NativeLinearEquationSolver.cpp
View File

@ -206,12 +206,19 @@ namespace storm {
return converged;
}
template<typename ValueType>
NativeLinearEquationSolver<ValueType>::WalkerChaeData::WalkerChaeData(storm::storage::SparseMatrix<ValueType> const& originalMatrix, std::vector<ValueType> const& originalB) : t(storm::utility::convertNumber<ValueType>(1000.0)) {
computeWalkerChaeMatrix(originalMatrix);
computeNewB(originalB);
precomputeAuxiliaryData();
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::computeWalkerChaeMatrix() const {
storm::storage::BitVector columnsWithNegativeEntries(this->A->getColumnCount());
void NativeLinearEquationSolver<ValueType>::WalkerChaeData::computeWalkerChaeMatrix(storm::storage::SparseMatrix<ValueType> const& originalMatrix) {
storm::storage::BitVector columnsWithNegativeEntries(originalMatrix.getColumnCount());
ValueType zero = storm::utility::zero<ValueType>();
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<ValueType> 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<ValueType>();
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<storm::storage::SparseMatrix<ValueType>>(builder.build());
matrix = builder.build();
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::WalkerChaeData::computeNewB(std::vector<ValueType> const& originalB) {
b = std::vector<ValueType>(originalB);
b.resize(matrix.getRowCount());
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::WalkerChaeData::precomputeAuxiliaryData() {
columnSums = std::vector<ValueType>(matrix.getColumnCount());
for (auto const& e : matrix) {
columnSums[e.getColumn()] += e.getValue();
}
newX.resize(matrix.getRowCount());
}
template<typename ValueType>
@ -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<WalkerChaeData>(*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<std::vector<ValueType>>(b);
walkerChaeB->resize(x.size());
}
// Choose a value for t in the algorithm.
ValueType t = storm::utility::convertNumber<ValueType>(1000);
// Precompute some data.
std::vector<ValueType> columnSums(x.size());
for (auto const& e : *walkerChaeMatrix) {
STORM_LOG_ASSERT(e.getValue() >= storm::utility::zero<ValueType>(), "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<ValueType> currentAx(x.size());
walkerChaeMatrix->multiplyWithVector(x, currentAx);
// Create an auxiliary vector that intermediately stores the result of the Walker-Chae step.
std::vector<ValueType> tmpX(x.size());
// Set up references to the x-vectors used in the iteration loop.
std::vector<ValueType>* currentX = &x;
std::vector<ValueType>* nextX = &tmpX;
std::vector<ValueType>* nextX = &walkerChaeData->newX;
std::vector<ValueType> 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<ValueType> 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<typename ValueType>
void NativeLinearEquationSolver<ValueType>::clearCache() const {
jacobiDecomposition.reset();
walkerChaeMatrix.reset();
walkerChaeData.reset();
LinearEquationSolver<ValueType>::clearCache();
}

20
src/storm/solver/NativeLinearEquationSolver.h
View File

@ -68,8 +68,6 @@ namespace storm {
virtual bool solveEquationsSOR(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType const& omega) const;
virtual bool solveEquationsJacobi(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
virtual bool solveEquationsWalkerChae(std::vector<ValueType>& x, std::vector<ValueType> 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::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>>> jacobiDecomposition;
mutable std::unique_ptr<storm::storage::SparseMatrix<ValueType>> walkerChaeMatrix;
mutable std::unique_ptr<std::vector<ValueType>> walkerChaeB;
struct WalkerChaeData {
WalkerChaeData(storm::storage::SparseMatrix<ValueType> const& originalMatrix, std::vector<ValueType> const& originalB);
void computeWalkerChaeMatrix(storm::storage::SparseMatrix<ValueType> const& originalMatrix);
void computeNewB(std::vector<ValueType> const& originalB);
void precomputeAuxiliaryData();
storm::storage::SparseMatrix<ValueType> matrix;
std::vector<ValueType> b;
ValueType t;
// Auxiliary data.
std::vector<ValueType> columnSums;
std::vector<ValueType> newX;
};
mutable std::unique_ptr<WalkerChaeData> walkerChaeData;
};
template<typename ValueType>

42
src/storm/storage/SparseMatrix.cpp
View File

@ -781,6 +781,18 @@ namespace storm {
}
}
template<typename ValueType>
std::vector<ValueType> SparseMatrix<ValueType>::getRowSumVector() const {
std::vector<ValueType> 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<typename ValueType>
ValueType SparseMatrix<ValueType>::getConstrainedRowSum(index_type row, storm::storage::BitVector const& constraint) const {
ValueType result = storm::utility::zero<ValueType>();
@ -1465,25 +1477,31 @@ namespace storm {
const_iterator it = this->begin();
const_iterator ite;
std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
typename std::vector<ValueType>::iterator resultIterator = result.begin();
typename std::vector<ValueType>::iterator resultIteratorEnd = result.end();
typename std::vector<ValueType>::const_iterator xIterator = x.begin();
typename std::vector<ValueType>::const_iterator columnSumsIterator = columnSums.begin();
for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++xIterator, ++columnSumsIterator) {
*resultIterator = storm::utility::zero<ValueType>();
// Clear all previous entries.
ValueType zero = storm::utility::zero<ValueType>();
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<Interval>::performWalkerChaeStep(std::vector<Interval> const& x, std::vector<Interval> const& columnSums, std::vector<Interval> const& b, std::vector<Interval> const& ax, std::vector<Interval>& result) const {
void SparseMatrix<Interval>::performWalkerChaeStep(std::vector<Interval> const& x, std::vector<Interval> const& rowSums, std::vector<Interval> const& b, std::vector<Interval> const& ax, std::vector<Interval>& result) const {
STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported.");
}
#endif

11
src/storm/storage/SparseMatrix.h
View File

@ -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<ValueType> 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<value_type> 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.

62
src/test/storm/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
View File

@ -194,38 +194,35 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) {
EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
}
// Does not converge any more. :(
// {
// matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
// matrixBuilder.addNextValue(0, 1, 1);
// matrixBuilder.addNextValue(1, 2, 1);
// matrixBuilder.addNextValue(2, 0, 1);
// storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
//
// storm::models::sparse::StateLabeling ap(3);
// ap.addLabel("a");
// ap.addLabelToState("a", 2);
//
// dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
//
// auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
// factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
// factory->getSettings().setOmega(0.9);
// storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
//
// std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
//
// std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
// storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
//
// EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
// EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
// EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
// }
{
matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
matrixBuilder.addNextValue(0, 1, 1);
matrixBuilder.addNextValue(1, 2, 1);
matrixBuilder.addNextValue(2, 0, 1);
storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
storm::models::sparse::StateLabeling ap(3);
ap.addLabel("a");
ap.addLabelToState("a", 2);
dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::WalkerChae);
storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
}
}
// Test is disabled as it does not converge any more. :(
TEST(DISABLED_NativeDtmcPrctlModelCheckerTest, LRA) {
TEST(NativeDtmcPrctlModelCheckerTest, LRA) {
storm::storage::SparseMatrixBuilder<double> matrixBuilder;
std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
@ -276,8 +273,9 @@ TEST(DISABLED_NativeDtmcPrctlModelCheckerTest, LRA) {
dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
factory->getSettings().setOmega(0.9);
factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::WalkerChae);
factory->getSettings().setPrecision(1e-7);
factory->getSettings().setMaximalNumberOfIterations(50000);
storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");

Write Preview
Loading…
Cancel
Save