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

Made GMRES and LSCG solution methods work for linear equation solving. Some further work on scheduler guessing. 

Former-commit-id: f6b538394a
tempestpy_adaptions
dehnert 12 years ago
parent
15542d46da
commit
d168b1848e
6 changed files with 144 additions and 23 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. 6
      examples/mdp/consensus/coin.pctl
  2. 1
      resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h
  3. 104
      src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
  4. 38
      src/solver/GmmxxLinearEquationSolver.h
  5. 14
      src/storage/SparseMatrix.h
  6. 4
      src/utility/Settings.cpp

6
examples/mdp/consensus/coin.pctl
View File

@ -1,11 +1,11 @@
// C1 (with probability 1, all N processes finish the protocol)
Pmin=? [ F finished ]
// Pmin=? [ F finished ]
// C2 (minimum probability that the protocol finishes with all coins equal to v) (v=1,2)
// Results are same for v=1 and v=2 by symmetry
// Analytic bound is (K-1)/(2*K)
Pmin=? [ F finished & all_coins_equal_0 ]
Pmin=? [ F finished & all_coins_equal_1 ]
//Pmin=? [ F finished & all_coins_equal_0 ]
//Pmin=? [ F finished & all_coins_equal_1 ]
// Max probability of finishing protocol with coins not all equal
Pmax=? [ F finished & !agree ]

1
resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h
View File

@ -74,6 +74,7 @@
#include "gmm_kernel.h"
#include "gmm_iter.h"
#include "gmm_modified_gram_schmidt.h"
#include "gmm_dense_Householder.h"
namespace gmm {

104
src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
View File

@ -8,6 +8,10 @@
#ifndef STORM_MODELCHECKER_PRCTL_SPARSEMDPPRCTLMODELCHECKER_H_
#define STORM_MODELCHECKER_PRCTL_SPARSEMDPPRCTLMODELCHECKER_H_
#include <vector>
#include <stack>
#include <fstream>
#include "src/modelchecker/prctl/AbstractModelChecker.h"
#include "src/solver/AbstractNondeterministicLinearEquationSolver.h"
#include "src/solver/GmmxxLinearEquationSolver.h"
@ -16,9 +20,6 @@
#include "src/utility/graph.h"
#include "src/utility/Settings.h"
#include <vector>
#include <stack>
namespace storm {
namespace modelchecker {
namespace prctl {
@ -287,6 +288,8 @@ namespace storm {
// Create resulting vector.
std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
std::vector<uint_fast64_t> guessedScheduler;
// Check whether we need to compute exact probabilities for some states.
if (this->getInitialStates().isDisjointFrom(maybeStates) || qualitative) {
if (qualitative) {
@ -313,7 +316,7 @@ namespace storm {
std::vector<Type> b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, submatrix.getRowCount());
// Create vector for results for maybe states.
std::vector<Type> x = this->getInitialValueIterationValues(minimize, submatrix, subNondeterministicChoiceIndices, b, statesWithProbability1, maybeStates);
std::vector<Type> x = this->getInitialValueIterationValues(minimize, submatrix, subNondeterministicChoiceIndices, b, statesWithProbability1, maybeStates, &guessedScheduler);
// Solve the corresponding system of equations.
if (linearEquationSolver != nullptr) {
@ -330,6 +333,44 @@ namespace storm {
storm::utility::vector::setVectorValues<Type>(*result, statesWithProbability0, storm::utility::constGetZero<Type>());
storm::utility::vector::setVectorValues<Type>(*result, statesWithProbability1, storm::utility::constGetOne<Type>());
// Output a scheduler for debug purposes.
std::vector<uint_fast64_t> myScheduler(this->getModel().getNumberOfStates());
this->computeTakenChoices(this->minimumOperatorStack.top(), false, *result, myScheduler, this->getModel().getNondeterministicChoiceIndices());
std::vector<uint_fast64_t> stateColoring(this->getModel().getNumberOfStates());
for (auto target : statesWithProbability1) {
stateColoring[target] = 1;
}
std::vector<std::string> colors(2);
colors[0] = "white";
colors[1] = "blue";
std::ofstream outFile;
outFile.open("real.dot");
storm::storage::BitVector filterStates(this->getModel().getNumberOfStates(), true);
this->getModel().writeDotToStream(outFile, true, &filterStates, result, nullptr, &stateColoring, &colors, &myScheduler);
outFile.close();
std::cout << "=========== Scheduler Comparison ===========" << std::endl;
uint_fast64_t index = 0;
for (auto state : maybeStates) {
if (myScheduler[state] != guessedScheduler[index]) {
std::cout << state << " right: " << myScheduler[state] << ", wrong: " << guessedScheduler[index] << std::endl;
std::cout << "Correct: " << std::endl;
typename storm::storage::SparseMatrix<Type>::Rows correctRow = this->getModel().getTransitionMatrix().getRow(this->getModel().getNondeterministicChoiceIndices()[state] + myScheduler[state]);
for (auto& transition : correctRow) {
std::cout << "target " << transition.column() << " with prob " << transition.value() << std::endl;
}
std::cout << "Incorrect: " << std::endl;
typename storm::storage::SparseMatrix<Type>::Rows incorrectRow = this->getModel().getTransitionMatrix().getRow(this->getModel().getNondeterministicChoiceIndices()[state] + myScheduler[state]);
for (auto& transition : incorrectRow) {
std::cout << "target " << transition.column() << " with prob " << transition.value() << std::endl;
}
}
++index;
}
// If we were required to generate a scheduler, do so now.
if (scheduler != nullptr) {
this->computeTakenChoices(this->minimumOperatorStack.top(), false, *result, *scheduler, this->getModel().getNondeterministicChoiceIndices());
@ -648,16 +689,36 @@ namespace storm {
std::vector<uint_fast64_t> const& subNondeterministicChoiceIndices,
std::vector<Type> const& rightHandSide,
storm::storage::BitVector const& targetStates,
storm::storage::BitVector const& maybeStates) const {
storm::storage::BitVector const& maybeStates,
std::vector<uint_fast64_t>* guessedScheduler = nullptr) const {
storm::settings::Settings* s = storm::settings::instance();
double precision = s->get<double>("precision");
if (s->get<bool>("use-heuristic-presolve")) {
std::pair<std::vector<Type>, std::vector<uint_fast64_t>> distancesAndPredecessorsPair = storm::utility::graph::performDijkstra(this->getModel(),
this->getModel().template getBackwardTransitions<Type>([](Type const& value) -> Type { return value; }),
minimize? ~(maybeStates | targetStates) : targetStates, &maybeStates);
minimize ? ~(maybeStates | targetStates) : targetStates, &maybeStates);
std::pair<std::vector<Type>, std::vector<uint_fast64_t>> distancesAndPredecessorsPair2 = storm::utility::graph::performDijkstra(this->getModel(),
this->getModel().template getBackwardTransitions<Type>([](Type const& value) -> Type { return value; }),
minimize ? targetStates : ~(maybeStates | targetStates), &maybeStates);
std::vector<uint_fast64_t> scheduler = this->convertShortestPathsToScheduler(false, maybeStates, distancesAndPredecessorsPair.first, distancesAndPredecessorsPair.second);
std::vector<uint_fast64_t> stateColoring(this->getModel().getNumberOfStates());
for (auto target : targetStates) {
stateColoring[target] = 1;
}
std::vector<std::string> colors(2);
colors[0] = "white";
colors[1] = "blue";
std::ofstream outFile;
outFile.open("guessed.dot");
storm::storage::BitVector filterStates(this->getModel().getNumberOfStates(), true);
this->getModel().writeDotToStream(outFile, true, &filterStates, &distancesAndPredecessorsPair.first, &distancesAndPredecessorsPair2.first, &stateColoring, &colors);
outFile.close();
std::vector<uint_fast64_t> scheduler = this->convertShortestPathsToScheduler(maybeStates, distancesAndPredecessorsPair.second);
std::vector<Type> result(scheduler.size(), Type(0.5));
std::vector<Type> b(scheduler.size());
storm::utility::vector::selectVectorValues(b, scheduler, subNondeterministicChoiceIndices, rightHandSide);
@ -674,28 +735,39 @@ namespace storm {
}
}
if (guessedScheduler != nullptr) {
*guessedScheduler = std::move(scheduler);
}
return result;
} else {
return std::vector<Type>(submatrix.getColumnCount());
}
}
std::vector<uint_fast64_t> convertShortestPathsToScheduler(storm::storage::BitVector const& maybeStates, std::vector<uint_fast64_t> const& shortestPathSuccessors) const {
std::vector<uint_fast64_t> convertShortestPathsToScheduler(bool minimize, storm::storage::BitVector const& maybeStates, std::vector<Type> const& distances, std::vector<uint_fast64_t> const& shortestPathSuccessors) const {
std::vector<uint_fast64_t> scheduler(maybeStates.getNumberOfSetBits());
Type maxProbability = 0;
Type extremalProbability = minimize ? 1 : 0;
Type currentProbability = 0;
uint_fast64_t currentStateIndex = 0;
for (auto state : maybeStates) {
maxProbability = 0;
extremalProbability = minimize ? 1 : 0;
for (uint_fast64_t row = 0, rowEnd = this->getModel().getNondeterministicChoiceIndices()[state + 1] - this->getModel().getNondeterministicChoiceIndices()[state]; row < rowEnd; ++row) {
typename storm::storage::SparseMatrix<Type>::Rows currentRow = this->getModel().getTransitionMatrix().getRows(this->getModel().getNondeterministicChoiceIndices()[state] + row, this->getModel().getNondeterministicChoiceIndices()[state] + row);
typename storm::storage::SparseMatrix<Type>::Rows currentRow = this->getModel().getTransitionMatrix().getRow(this->getModel().getNondeterministicChoiceIndices()[state] + row);
currentProbability = 0;
for (auto& transition : currentRow) {
if (transition.column() == shortestPathSuccessors[state] && transition.value() > maxProbability) {
maxProbability = transition.value();
scheduler[currentStateIndex] = row;
}
currentProbability += transition.value() * (1 / std::exp(distances[transition.column()]));
}
if (minimize && currentProbability < extremalProbability) {
extremalProbability = currentProbability;
scheduler[currentStateIndex] = row;
} else if (!minimize && currentProbability > extremalProbability) {
extremalProbability = currentProbability;
scheduler[currentStateIndex] = row;
}
}

38
src/solver/GmmxxLinearEquationSolver.h
View File

@ -86,6 +86,44 @@ namespace storm {
gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter);
}
// Check if the solver converged and issue a warning otherwise.
if (iter.converged()) {
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
} else {
LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
}
delete gmmA;
} else if (s->getString("lemethod") == "lscg") {
LOG4CPLUS_INFO(logger, "Using LSCG method.");
// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
if (precond != "none") {
LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the LSCG method. Dropping preconditioner.");
}
gmm::least_squares_cg(*gmmA, x, b, iter);
// Check if the solver converged and issue a warning otherwise.
if (iter.converged()) {
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
} else {
LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
}
delete gmmA;
} else if (s->getString("lemethod") == "gmres") {
LOG4CPLUS_INFO(logger, "Using GMRES method.");
// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
if (precond == "ilu") {
gmm::gmres(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
} else if (precond == "diagonal") {
gmm::gmres(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
} else if (precond == "ildlt") {
gmm::gmres(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
} else if (precond == "none") {
gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), 50, iter);
}
// Check if the solver converged and issue a warning otherwise.
if (iter.converged()) {
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");

14
src/storage/SparseMatrix.h
View File

@ -1170,14 +1170,24 @@ public:
/*!
* Returns an object representing the consecutive rows given by the parameters.
*
* @param The starting row.
* @param The ending row (which is included in the result).
* @param startRow The starting row.
* @param endRow The ending row (which is included in the result).
* @return An object representing the consecutive rows given by the parameters.
*/
Rows getRows(uint_fast64_t startRow, uint_fast64_t endRow) const {
return Rows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]);
}
/*!
* Returns an object representing the given row.
*
* @param row The chosen row.
* @return An object representing the given row.
*/
Rows getRow(uint_fast64_t row) const {
return getRows(row, row);
}
/*!
* Returns a const iterator to the rows of the matrix.
*

4
src/utility/Settings.cpp
View File

@ -129,7 +129,7 @@ void checkExplicit(const std::vector<std::string>& filenames) {
* Throws an exception of type InvalidSettings in case the selected method is illegal.
*/
static void validateLeMethod(const std::string& lemethod) {
if ((lemethod != "bicgstab") && (lemethod != "qmr") && (lemethod != "jacobi")) {
if ((lemethod != "bicgstab") && (lemethod != "qmr") && (lemethod != "jacobi") && (lemethod != "lscg") && (lemethod != "gmres")) {
throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid.";
}
}
@ -165,7 +165,7 @@ void Settings::initDescriptions() {
("transrew", bpo::value<std::string>()->default_value(""), "name of transition reward file")
("staterew", bpo::value<std::string>()->default_value(""), "name of state reward file")
("fix-deadlocks", "insert self-loops for states without outgoing transitions")
("lemethod", boost::program_options::value<std::string>()->default_value("bicgstab")->notifier(&storm::settings::validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}.")
("lemethod", boost::program_options::value<std::string>()->default_value("bicgstab")->notifier(&storm::settings::validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, lscg, gmres, jacobi}.")
("maxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving.")
("precision", boost::program_options::value<double>()->default_value(1e-6), "Sets the precision for iterative equation solving.")
("precond", boost::program_options::value<std::string>()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}.")

Write Preview
Loading…
Cancel
Save