diff --git a/examples/mdp/consensus/coin.pctl b/examples/mdp/consensus/coin.pctl index 7eb1c2bc8..224ffc495 100644 --- a/examples/mdp/consensus/coin.pctl +++ b/examples/mdp/consensus/coin.pctl @@ -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 ] diff --git a/resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h b/resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h index 38f760ae4..caaa37e2a 100644 --- a/resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h +++ b/resources/3rdparty/gmm-4.2/include/gmm/gmm_solver_gmres.h @@ -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 { diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 66a0fb794..3c04bcfbe 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -8,6 +8,10 @@ #ifndef STORM_MODELCHECKER_PRCTL_SPARSEMDPPRCTLMODELCHECKER_H_ #define STORM_MODELCHECKER_PRCTL_SPARSEMDPPRCTLMODELCHECKER_H_ +#include +#include +#include + #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 -#include - namespace storm { namespace modelchecker { namespace prctl { @@ -287,6 +288,8 @@ namespace storm { // Create resulting vector. std::vector* result = new std::vector(this->getModel().getNumberOfStates()); + std::vector 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 b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, submatrix.getRowCount()); // Create vector for results for maybe states. - std::vector x = this->getInitialValueIterationValues(minimize, submatrix, subNondeterministicChoiceIndices, b, statesWithProbability1, maybeStates); + std::vector 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(*result, statesWithProbability0, storm::utility::constGetZero()); storm::utility::vector::setVectorValues(*result, statesWithProbability1, storm::utility::constGetOne()); + // Output a scheduler for debug purposes. + std::vector myScheduler(this->getModel().getNumberOfStates()); + this->computeTakenChoices(this->minimumOperatorStack.top(), false, *result, myScheduler, this->getModel().getNondeterministicChoiceIndices()); + + std::vector stateColoring(this->getModel().getNumberOfStates()); + for (auto target : statesWithProbability1) { + stateColoring[target] = 1; + } + + std::vector 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::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::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 const& subNondeterministicChoiceIndices, std::vector const& rightHandSide, storm::storage::BitVector const& targetStates, - storm::storage::BitVector const& maybeStates) const { + storm::storage::BitVector const& maybeStates, + std::vector* guessedScheduler = nullptr) const { storm::settings::Settings* s = storm::settings::instance(); double precision = s->get("precision"); if (s->get("use-heuristic-presolve")) { - std::pair, std::vector> distancesAndPredecessorsPair = storm::utility::graph::performDijkstra(this->getModel(), this->getModel().template getBackwardTransitions([](Type const& value) -> Type { return value; }), - minimize? ~(maybeStates | targetStates) : targetStates, &maybeStates); + minimize ? ~(maybeStates | targetStates) : targetStates, &maybeStates); + + std::pair, std::vector> distancesAndPredecessorsPair2 = storm::utility::graph::performDijkstra(this->getModel(), + this->getModel().template getBackwardTransitions([](Type const& value) -> Type { return value; }), + minimize ? targetStates : ~(maybeStates | targetStates), &maybeStates); + + std::vector scheduler = this->convertShortestPathsToScheduler(false, maybeStates, distancesAndPredecessorsPair.first, distancesAndPredecessorsPair.second); + + std::vector stateColoring(this->getModel().getNumberOfStates()); + for (auto target : targetStates) { + stateColoring[target] = 1; + } + + std::vector 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 scheduler = this->convertShortestPathsToScheduler(maybeStates, distancesAndPredecessorsPair.second); std::vector result(scheduler.size(), Type(0.5)); std::vector 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(submatrix.getColumnCount()); } } - std::vector convertShortestPathsToScheduler(storm::storage::BitVector const& maybeStates, std::vector const& shortestPathSuccessors) const { + std::vector convertShortestPathsToScheduler(bool minimize, storm::storage::BitVector const& maybeStates, std::vector const& distances, std::vector const& shortestPathSuccessors) const { std::vector 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::Rows currentRow = this->getModel().getTransitionMatrix().getRows(this->getModel().getNondeterministicChoiceIndices()[state] + row, this->getModel().getNondeterministicChoiceIndices()[state] + row); + typename storm::storage::SparseMatrix::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; } } diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 2f601f3ed..306ba38b8 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -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* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(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* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + if (precond == "ilu") { + gmm::gmres(*gmmA, x, b, gmm::ilu_precond>(*gmmA), 50, iter); + } else if (precond == "diagonal") { + gmm::gmres(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), 50, iter); + } else if (precond == "ildlt") { + gmm::gmres(*gmmA, x, b, gmm::ildlt_precond>(*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."); diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index d57fd2e23..f0bb0e0f1 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -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. * diff --git a/src/utility/Settings.cpp b/src/utility/Settings.cpp index 6fbfa07de..97e50e562 100644 --- a/src/utility/Settings.cpp +++ b/src/utility/Settings.cpp @@ -129,7 +129,7 @@ void checkExplicit(const std::vector& 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()->default_value(""), "name of transition reward file") ("staterew", bpo::value()->default_value(""), "name of state reward file") ("fix-deadlocks", "insert self-loops for states without outgoing transitions") - ("lemethod", boost::program_options::value()->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()->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()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving.") ("precision", boost::program_options::value()->default_value(1e-6), "Sets the precision for iterative equation solving.") ("precond", boost::program_options::value()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}.")