From abf6f85b63ca7db2a3f749306728c92741b73591 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 18 Jun 2013 17:30:49 +0200 Subject: [PATCH] Intermediate commit to switch workplace. Former-commit-id: 11932e19d7a727ff181bb2dc0389ab370adad87e --- CMakeLists.txt | 2 + .../prctl/GmmxxDtmcPrctlModelChecker.h | 309 ---- .../prctl/GmmxxMdpPrctlModelChecker.h | 186 --- .../prctl/SparseDtmcPrctlModelChecker.h | 71 +- .../prctl/SparseMdpPrctlModelChecker.h | 1315 ++++++++--------- src/models/AtomicPropositionsLabeling.h | 17 + src/solver/AbstractLinearEquationSolver.h | 26 + ...ractNondeterministicLinearEquationSolver.h | 175 +++ src/solver/GmmxxLinearEquationSolver.h | 228 +++ ...mmxxNondeterministicLinearEquationSolver.h | 122 ++ src/storage/BitVector.h | 20 + src/storm.cpp | 9 +- src/utility/Settings.cpp | 10 +- src/utility/Settings.h | 10 +- src/utility/vector.h | 17 + 15 files changed, 1270 insertions(+), 1247 deletions(-) delete mode 100644 src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h delete mode 100644 src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h create mode 100644 src/solver/AbstractLinearEquationSolver.h create mode 100644 src/solver/AbstractNondeterministicLinearEquationSolver.h create mode 100644 src/solver/GmmxxLinearEquationSolver.h create mode 100644 src/solver/GmmxxNondeterministicLinearEquationSolver.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 874f2133c..a85cc4757 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -141,6 +141,7 @@ file(GLOB_RECURSE STORM_STORAGE_FILES ${PROJECT_SOURCE_DIR}/src/storage/*.h ${PR file(GLOB_RECURSE STORM_UTILITY_FILES ${PROJECT_SOURCE_DIR}/src/utility/*.h ${PROJECT_SOURCE_DIR}/src/utility/*.cpp) file(GLOB STORM_IR_FILES ${PROJECT_SOURCE_DIR}/src/ir/*.h ${PROJECT_SOURCE_DIR}/src/ir/*.cpp) file(GLOB_RECURSE STORM_IR_EXPRESSIONS_FILES ${PROJECT_SOURCE_DIR}/src/ir/expressions/*.h ${PROJECT_SOURCE_DIR}/src/ir/expressions/*.cpp) +file(GLOB_RECURSE STORM_SOLVER_FILES ${PROJECT_SOURCE_DIR}/src/solver/*.h ${PROJECT_SOURCE_DIR}/src/solver/*.cpp) # Test Sources # Note that the tests also need the source files, except for the main file @@ -164,6 +165,7 @@ source_group(storage FILES ${STORM_STORAGE_FILES}) source_group(utility FILES ${STORM_UTILITY_FILES}) source_group(ir FILES ${STORM_IR_FILES}) source_group(ir\\expressions FILES ${STORM_IR_EXPRESSIONS_FILES}) +source_group(solver FILES ${STORM_SOLVER_FILES}) source_group(functional-test FILES ${STORM_FUNCTIONAL_TEST_FILES}) source_group(performance-test FILES ${STORM_PERFORMANCE_TEST_FILES}) diff --git a/src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h b/src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h deleted file mode 100644 index fbbaed0ea..000000000 --- a/src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h +++ /dev/null @@ -1,309 +0,0 @@ -/* - * GmmxxDtmcPrctlModelChecker.h - * - * Created on: 06.12.2012 - * Author: Christian Dehnert - */ - -#ifndef STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ -#define STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ - -#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" -#include "src/adapters/GmmxxAdapter.h" -#include "src/storage/SparseMatrix.h" -#include "src/utility/ConstTemplates.h" -#include "src/utility/Settings.h" - -#include "gmm/gmm_matrix.h" -#include "gmm/gmm_iter_solvers.h" - -#include - -namespace storm { -namespace modelchecker { -namespace prctl { - -/* - * An implementation of the SparseDtmcPrctlModelChecker interface that uses gmm++ as the solving backend. - */ -template -class GmmxxDtmcPrctlModelChecker : public SparseDtmcPrctlModelChecker { - -public: - /*! - * Constructs a GmmxxDtmcPrctlModelChecker with the given model. - * - * @param model The DTMC to be checked. - */ - explicit GmmxxDtmcPrctlModelChecker(storm::models::Dtmc const& dtmc) : SparseDtmcPrctlModelChecker(dtmc) { - // Intentionally left empty. - } - - /*! - * Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly - * constructed model checker will have the model of the given model checker as its associated model. - */ - explicit GmmxxDtmcPrctlModelChecker(storm::modelchecker::prctl::GmmxxDtmcPrctlModelChecker const& modelchecker) : SparseDtmcPrctlModelChecker(modelchecker) { - // Intentionally left empty. - } - - /*! - * Virtual destructor. Needs to be virtual, because this class has virtual methods. - */ - virtual ~GmmxxDtmcPrctlModelChecker() { - // Intentionally left empty. - } - - /*! - * Returns the name of this module. - * @returns The name of this module. - */ - static std::string getModuleName() { - return "gmm++det"; - } - - /*! - * Returns a trigger such that if the option "matrixlib" is set to "gmm++", this model checker - * is to be used. - * @returns An option trigger for this module. - */ - static std::pair getOptionTrigger() { - return std::pair("matrixlib", "gmm++"); - } - - /*! - * Registers all options associated with the gmm++ matrix library. - */ - static void putOptions(boost::program_options::options_description* desc) { - desc->add_options()("lemethod", boost::program_options::value()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}."); - desc->add_options()("maxiter", boost::program_options::value()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); - desc->add_options()("precision", boost::program_options::value()->default_value(1e-6), "Sets the precision for iterative equation solving."); - desc->add_options()("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}."); - desc->add_options()("relative", boost::program_options::value()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence."); - } - - /*! - * Validates whether the given lemethod matches one of the available ones. - * 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")) { - throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid."; - } - } - - /*! - * Validates whether the given preconditioner matches one of the available ones. - * Throws an exception of type InvalidSettings in case the selected preconditioner is illegal. - */ - static void validatePreconditioner(const std::string& preconditioner) { - if ((preconditioner != "ilu") && (preconditioner != "diagonal") && (preconditioner != "ildlt") && (preconditioner != "none")) { - throw exceptions::InvalidSettingsException() << "Argument " << preconditioner << " for option 'precond' is invalid."; - } - } - -private: - - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b, uint_fast64_t n = 1) const { - // Transform the transition probability A to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - - // Set up some temporary variables so that we can just swap pointers instead of copying the result after each - // iteration. - std::vector* swap = nullptr; - std::vector* currentVector = &x; - std::vector* tmpVector = new std::vector(this->getModel().getNumberOfStates()); - - // Now perform matrix-vector multiplication as long as we meet the bound. - for (uint_fast64_t i = 0; i < n; ++i) { - gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector); - swap = tmpVector; - tmpVector = currentVector; - currentVector = swap; - - // If requested, add an offset to the current result vector. - if (b != nullptr) { - gmm::add(*b, *currentVector); - } - } - - // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, because - // the output is supposed to be stored in x. - if (n % 2 == 1) { - std::swap(x, *currentVector); - delete currentVector; - } else { - delete tmpVector; - } - - delete gmmxxMatrix; - } - - virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - // Get the settings object to customize linear solving. - storm::settings::Settings* s = storm::settings::instance(); - - // Prepare an iteration object that determines the accuracy, maximum number of iterations - // and the like. - gmm::iteration iter(s->get("precision"), 0, s->get("maxiter")); - - // Print some information about the used preconditioner. - const std::string& precond = s->getString("precond"); - LOG4CPLUS_INFO(logger, "Starting iterative solver."); - if (s->getString("lemethod") == "jacobi") { - if (precond != "none") { - LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the Jacobi method. Dropping preconditioner."); - } - } else { - if (precond == "ilu") { - LOG4CPLUS_INFO(logger, "Using ILU preconditioner."); - } else if (precond == "diagonal") { - LOG4CPLUS_INFO(logger, "Using diagonal preconditioner."); - } else if (precond == "ildlt") { - LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner."); - } else if (precond == "none") { - LOG4CPLUS_INFO(logger, "Using no preconditioner."); - } - } - - // Now do the actual solving. - if (s->getString("lemethod") == "bicgstab") { - LOG4CPLUS_INFO(logger, "Using BiCGStab 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::bicgstab(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); - } else if (precond == "diagonal") { - gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); - } else if (precond == "ildlt") { - gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); - } else if (precond == "none") { - gmm::bicgstab(*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") == "qmr") { - LOG4CPLUS_INFO(logger, "Using QMR 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::qmr(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); - } else if (precond == "diagonal") { - gmm::qmr(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); - } else if (precond == "ildlt") { - gmm::qmr(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); - } else if (precond == "none") { - 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") == "jacobi") { - LOG4CPLUS_INFO(logger, "Using Jacobi method."); - uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b); - - // Check if the solver converged and issue a warning otherwise. - if (iterations < s->get("maxiter")) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - } - } - - /*! - * Solves the linear equation system A*x = b given by the parameters using the Jacobi method. - * - * @param A The matrix specifying the coefficients of the linear equations. - * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but - * may be ignored. - * @param b The right-hand side of the equation system. - * @returns The solution vector x of the system of linear equations as the content of the parameter x. - * @returns The number of iterations needed until convergence. - */ - uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - // Get the settings object to customize linear solving. - storm::settings::Settings* s = storm::settings::instance(); - - double precision = s->get("precision"); - uint_fast64_t maxIterations = s->get("maxiter"); - bool relative = s->get("relative"); - - // Get a Jacobi decomposition of the matrix A. - storm::storage::SparseMatrix::SparseJacobiDecomposition_t jacobiDecomposition = A.getJacobiDecomposition(); - - // Convert the (inverted) diagonal matrix to gmm++'s format. - gmm::csr_matrix* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.second)); - // Convert the LU matrix to gmm++'s format. - gmm::csr_matrix* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); - - LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver."); - - // x_(k + 1) = D^-1 * (b - R * x_k) - // In x we keep a copy of the result for swapping in the loop (e.g. less copy-back). - std::vector* xNext = new std::vector(x.size()); - const std::vector* xCopy = xNext; - std::vector* xCurrent = &x; - - // Target vector for precision calculation. - std::vector* residuum = new std::vector(x.size()); - - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; - - while (!converged && iterationCount < maxIterations) { - // R * x_k (xCurrent is x_k) -> xNext - gmm::mult(*gmmxxLU, *xCurrent, *xNext); - // b - R * x_k (xNext contains R * x_k) -> xNext - gmm::add(b, gmm::scaled(*xNext, -1.0), *xNext); - // D^-1 * (b - R * x_k) -> xNext - gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext); - - // Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the - // vector. - std::vector *const swap = xNext; - xNext = xCurrent; - xCurrent = swap; - - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*xCurrent, *xNext, precision, relative); - - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } - - // 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 parameter x. - if (xCurrent == xCopy) { - std::swap(x, *xCurrent); - } - - // As xCopy always points to the copy of x used for swapping, we can safely delete it. - delete xCopy; - - // Also delete the other dynamically allocated variables. - delete residuum; - delete gmmxxDiagonalInverted; - delete gmmxxLU; - - return iterationCount; - } -}; - -} // namespace prctl -} // namespace modelchecker -} // namespace storm - -#endif /* STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ */ diff --git a/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h deleted file mode 100644 index ff3ebd971..000000000 --- a/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h +++ /dev/null @@ -1,186 +0,0 @@ -/* - * GmmxxDtmcPrctlModelChecker.h - * - * Created on: 06.12.2012 - * Author: Christian Dehnert - */ - -#ifndef STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ -#define STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ - -#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" -#include "src/adapters/GmmxxAdapter.h" - -#include "gmm/gmm_matrix.h" -#include "gmm/gmm_iter_solvers.h" - -#include - -namespace storm { -namespace modelchecker { -namespace prctl { - -/* - * An implementation of the SparseMdpPrctlModelChecker interface that uses gmm++ as the solving backend. - */ -template -class GmmxxMdpPrctlModelChecker : public SparseMdpPrctlModelChecker { - -public: - /*! - * Constructs a GmmxxMdpPrctlModelChecker with the given model. - * - * @param model The MDP to be checked. - */ - explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp const& model) : SparseMdpPrctlModelChecker(model) { - // Intentionally left empty. - } - - /*! - * Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly - * constructed model checker will have the model of the given model checker as its associated model. - */ - explicit GmmxxMdpPrctlModelChecker(storm::modelchecker::prctl::GmmxxMdpPrctlModelChecker const& modelchecker) : SparseMdpPrctlModelChecker(modelchecker) { - // Intentionally left empty. - } - - /*! - * Virtual destructor. Needs to be virtual, because this class has virtual methods. - */ - virtual ~GmmxxMdpPrctlModelChecker() { - // Intentionally left empty. - } - -private: - - /*! - * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b - * until x[n], where x[0] = x. - * - * @param A The matrix that is to be multiplied against the vector. - * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, - * i.e. after the method returns, this vector will contain the computed values. - * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. - * @param b If not null, this vector is being added to the result after each matrix-vector multiplication. - * @param n Specifies the number of iterations the matrix-vector multiplication is performed. - * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. - */ - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& nondeterministicChoiceIndices, std::vector* b = nullptr, uint_fast64_t n = 1) const { - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - - // Create vector for result of multiplication, which is reduced to the result vector after - // each multiplication. - std::vector multiplyResult(A.getRowCount()); - - // Now perform matrix-vector multiplication as long as we meet the bound of the formula. - for (uint_fast64_t i = 0; i < n; ++i) { - gmm::mult(*gmmxxMatrix, x, multiplyResult); - - if (b != nullptr) { - gmm::add(*b, multiplyResult); - } - - if (this->minimumOperatorStack.top()) { - storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); - } else { - storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); - } - } - - // Delete intermediate results and return result. - delete gmmxxMatrix; - } - - /*! - * Solves the given equation system under the given parameters using the power method. - * - * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. - * @param A The matrix A specifying the coefficients of the equations. - * @param x The vector x for which to solve the equations. The initial value of the elements of - * this vector are used as the initial guess and might thus influence performance and convergence. - * @param b The vector b specifying the values on the right-hand-sides of the equations. - * @return The solution of the system of linear equations in form of the elements of the vector - * x. - */ - void solveEquationSystem(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector const& nondeterministicChoiceIndices, std::vector* takenChoices = nullptr) const override { - // Get the settings object to customize solving. - storm::settings::Settings* s = storm::settings::instance(); - - // Get relevant user-defined settings for solving the equations. - double precision = s->get("precision"); - unsigned maxIterations = s->get("maxiter"); - bool relative = s->get("relative"); - - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - - // Set up the environment for the power method. - std::vector multiplyResult(A.getRowCount()); - std::vector* currentX = &x; - std::vector* newX = new std::vector(x.size()); - std::vector* swap = nullptr; - uint_fast64_t iterations = 0; - bool converged = false; - - // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number - // of iterations. - while (!converged && iterations < maxIterations) { - std::cout << "iter: " << iterations << std::endl; - // Compute x' = A*x + b. - gmm::mult(*gmmxxMatrix, *currentX, multiplyResult); - gmm::add(b, multiplyResult); - - // Reduce the vector x' by applying min/max for all non-deterministic choices. - if (minimize) { - storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); - } else { - storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); - } - - // Determine whether the method converged. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); - - if (!converged) { - for (uint_fast64_t i = 0; i < newX->size(); ++i) { - std::cout << (*currentX)[i] << " vs. " << (*newX)[i] << std::endl; - } - } - - - // Update environment variables. - swap = currentX; - currentX = newX; - newX = swap; - ++iterations; - } - - // If we were requested to record the taken choices, we have to construct the vector now. - if (takenChoices != nullptr) { - this->computeTakenChoices(minimize, multiplyResult, *takenChoices, nondeterministicChoiceIndices); - } - - // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result - // is currently stored in currentX, but x is the output vector. - if (iterations % 2 == 1) { - std::swap(x, *currentX); - delete currentX; - } else { - delete newX; - } - delete gmmxxMatrix; - - // Check if the solver converged and issue a warning otherwise. - if (converged) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); - } - } -}; - -} //namespace prctl -} //namespace modelchecker -} //namespace storm - -#endif /* STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ */ diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index 7c70e59df..50ab036c1 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -10,6 +10,7 @@ #include "src/modelchecker/prctl/AbstractModelChecker.h" #include "src/models/Dtmc.h" +#include "src/solver/AbstractLinearEquationSolver.h" #include "src/utility/vector.h" #include "src/utility/graph.h" @@ -26,21 +27,22 @@ namespace prctl { template class SparseDtmcPrctlModelChecker : public AbstractModelChecker { -public: /*! * Constructs a SparseDtmcPrctlModelChecker with the given model. * * @param model The DTMC to be checked. */ - explicit SparseDtmcPrctlModelChecker(storm::models::Dtmc const& model) : AbstractModelChecker(model) { + explicit SparseDtmcPrctlModelChecker(storm::models::Dtmc const& model, storm::solver::AbstractLinearEquationSolver* linearEquationSolver) : AbstractModelChecker(model), linearEquationSolver(linearEquationSolver) { // Intentionally left empty. } + +public: /*! * Copy constructs a SparseDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly * constructed model checker will have the model of the given model checker as its associated model. */ - explicit SparseDtmcPrctlModelChecker(storm::modelchecker::prctl::SparseDtmcPrctlModelChecker const& modelChecker) : AbstractModelChecker(modelChecker) { + explicit SparseDtmcPrctlModelChecker(storm::modelchecker::prctl::SparseDtmcPrctlModelChecker const& modelChecker) : AbstractModelChecker(modelChecker), linearEquationSolver(modelChecker.linearEquationSolver->clone()) { // Intentionally left empty. } @@ -120,7 +122,11 @@ public: storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne()); // Perform the matrix vector multiplication as often as required by the formula bound. - this->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound()); + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } // Set the values of the resulting vector accordingly. storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, subresult); @@ -156,8 +162,12 @@ public: delete nextStates; // Perform one single matrix-vector multiplication. - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result); - + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + // Return result. return result; } @@ -282,7 +292,11 @@ public: std::vector b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); // Now solve the created system of linear equations. - this->solveEquationSystem(submatrix, x, b); + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->solveEquationSystem(submatrix, x, b); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } // Set values of resulting vector according to result. storm::utility::vector::setVectorValues(*result, maybeStates, x); @@ -317,7 +331,11 @@ public: std::vector* result = new std::vector(this->getModel().getStateRewardVector()); // Perform the actual matrix-vector multiplication as long as the bound of the formula is met. - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound()); + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } // Return result. return result; @@ -361,7 +379,11 @@ public: } // Perform the actual matrix-vector multiplication as long as the bound of the formula is met. - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound()); + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } // Return result. return result; @@ -448,7 +470,11 @@ public: } // Now solve the resulting equation system. - this->solveEquationSystem(submatrix, x, b); + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->solveEquationSystem(submatrix, x, b); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } // Set values of resulting vector according to result. storm::utility::vector::setVectorValues(*result, maybeStates, x); @@ -464,29 +490,8 @@ public: } private: - /*! - * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b - * until x[n], where x[0] = x. - * - * @param A The matrix that is to be multiplied against the vector. - * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, - * i.e. after the method returns, this vector will contain the computed values. - * @param b If not null, this vector is being added to the result after each matrix-vector multiplication. - * @param n Specifies the number of iterations the matrix-vector multiplication is performed. - * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. - */ - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1) const = 0; - - /*! - * Solves the equation system A*x = b given by the parameters. - * - * @param A The matrix specifying the coefficients of the linear equations. - * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but - * may be ignored. - * @param b The right-hand side of the equation system. - * @returns The solution vector x of the system of linear equations as the content of the parameter x. - */ - virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const = 0; + // An object that is used for solving linear equations and performing matrix-vector multiplication. + std::unique_ptr> linearEquationSolver; }; } // namespace prctl diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index eea889769..42173d0b7 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -10,6 +10,7 @@ #include "src/modelchecker/prctl/AbstractModelChecker.h" #include "src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h" +#include "src/solver/AbstractNondeterministicLinearEquationSolver.h" #include "src/models/Mdp.h" #include "src/utility/vector.h" #include "src/utility/graph.h" @@ -18,711 +19,619 @@ #include namespace storm { -namespace modelchecker { -namespace prctl { - -/*! - * @brief - * Interface for all model checkers that can verify PRCTL formulae over MDPs represented as a sparse matrix. - */ -template -class SparseMdpPrctlModelChecker : public AbstractModelChecker { - -public: - /*! - * Constructs a SparseMdpPrctlModelChecker with the given model. - * - * @param model The MDP to be checked. - */ - explicit SparseMdpPrctlModelChecker(storm::models::Mdp const& model) : AbstractModelChecker(model), minimumOperatorStack() { - // Intentionally left empty. - } - - /*! - * Copy constructs a SparseMdpPrctlModelChecker from the given model checker. In particular, this means that the newly - * constructed model checker will have the model of the given model checker as its associated model. - */ - explicit SparseMdpPrctlModelChecker(storm::modelchecker::prctl::SparseMdpPrctlModelChecker const& modelchecker) - : AbstractModelChecker(modelchecker), minimumOperatorStack() { - // Intentionally left empty. - } - - /*! - * Virtual destructor. Needs to be virtual, because this class has virtual methods. - */ - virtual ~SparseMdpPrctlModelChecker() { - // Intentionally left empty. - } - - /*! - * Returns a constant reference to the MDP associated with this model checker. - * @returns A constant reference to the MDP associated with this model checker. - */ - storm::models::Mdp const& getModel() const { - return AbstractModelChecker::template getModel>(); - } - - /*! - * Checks the given formula that is a P/R operator without a bound. - * - * @param formula The formula to check. - * @returns The set of states satisfying the formula represented by a bit vector. - */ - virtual std::vector* checkNoBoundOperator(const storm::property::prctl::AbstractNoBoundOperator& formula) const { - // Check if the operator was an non-optimality operator and report an error in that case. - if (!formula.isOptimalityOperator()) { - LOG4CPLUS_ERROR(logger, "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models."); - throw storm::exceptions::InvalidArgumentException() << "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models."; - } - minimumOperatorStack.push(formula.isMinimumOperator()); - std::vector* result = formula.check(*this, false); - minimumOperatorStack.pop(); - return result; - } - - /*! - * Checks the given formula that is a bounded-until formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @returns The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkBoundedUntil(const storm::property::prctl::BoundedUntil& formula, bool qualitative) const { - // First, we need to compute the states that satisfy the sub-formulas of the until-formula. - storm::storage::BitVector* leftStates = formula.getLeft().check(*this); - storm::storage::BitVector* rightStates = formula.getRight().check(*this); - std::vector* result = new std::vector(this->getModel().getNumberOfStates()); - - // Determine the states that have 0 probability of reaching the target states. - storm::storage::BitVector statesWithProbabilityGreater0; - if (this->minimumOperatorStack.top()) { - statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound()); - } else { - statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound()); - } - - // Check if we already know the result (i.e. probability 0) for all initial states and - // don't compute anything in this case. - if (this->getInitialStates().isDisjointFrom(statesWithProbabilityGreater0)) { - LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step." - << " No exact probabilities were computed."); - // Set the values for all maybe-states to 0.5 to indicate that their probability values are not 0 (and - // not necessarily 1). - storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, Type(0.5)); - } else { - // In this case we have have to compute the probabilities. - - // We can eliminate the rows and columns from the original transition probability matrix that have probability 0. - storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0, this->getModel().getNondeterministicChoiceIndices()); - - // Get the "new" nondeterministic choice indices for the submatrix. - std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(statesWithProbabilityGreater0); - - // Compute the new set of target states in the reduced system. - storm::storage::BitVector rightStatesInReducedSystem = statesWithProbabilityGreater0 % *rightStates; + namespace modelchecker { + namespace prctl { - // Make all rows absorbing that satisfy the second sub-formula. - submatrix.makeRowsAbsorbing(rightStatesInReducedSystem, subNondeterministicChoiceIndices); + /*! + * This class represents the base class for all PRCTL model checkers for MDPs. + */ + template + class SparseMdpPrctlModelChecker : public AbstractModelChecker { + + public: + /*! + * Constructs a SparseMdpPrctlModelChecker with the given model. + * + * @param model The MDP to be checked. + */ + explicit SparseMdpPrctlModelChecker(storm::models::Mdp const& model, storm::solver::AbstractNondeterministicLinearEquationSolver* linearEquationSolver) : AbstractModelChecker(model), minimumOperatorStack(), linearEquationSolver(linearEquationSolver) { + // Intentionally left empty. + } + + /*! + * Copy constructs a SparseMdpPrctlModelChecker from the given model checker. In particular, this means that the newly + * constructed model checker will have the model of the given model checker as its associated model. + */ + explicit SparseMdpPrctlModelChecker(storm::modelchecker::prctl::SparseMdpPrctlModelChecker const& modelchecker) + : AbstractModelChecker(modelchecker), minimumOperatorStack(), linearEquationSolver(new storm::solver::AbstractNondeterministicLinearEquationSolver()) { + // Intentionally left empty. + } + + /*! + * Returns a constant reference to the MDP associated with this model checker. + * @returns A constant reference to the MDP associated with this model checker. + */ + storm::models::Mdp const& getModel() const { + return AbstractModelChecker::template getModel>(); + } + + /*! + * Checks the given formula that is a P/R operator without a bound. + * + * @param formula The formula to check. + * @returns The set of states satisfying the formula represented by a bit vector. + */ + virtual std::vector* checkNoBoundOperator(const storm::property::prctl::AbstractNoBoundOperator& formula) const { + // Check if the operator was an non-optimality operator and report an error in that case. + if (!formula.isOptimalityOperator()) { + LOG4CPLUS_ERROR(logger, "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models."); + throw storm::exceptions::InvalidArgumentException() << "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models."; + } + minimumOperatorStack.push(formula.isMinimumOperator()); + std::vector* result = formula.check(*this, false); + minimumOperatorStack.pop(); + return result; + } + + /*! + * Checks the given formula that is a bounded-until formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @returns The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkBoundedUntil(const storm::property::prctl::BoundedUntil& formula, bool qualitative) const { + // First, we need to compute the states that satisfy the sub-formulas of the until-formula. + storm::storage::BitVector* leftStates = formula.getLeft().check(*this); + storm::storage::BitVector* rightStates = formula.getRight().check(*this); + std::vector* result = new std::vector(this->getModel().getNumberOfStates()); + + // Determine the states that have 0 probability of reaching the target states. + storm::storage::BitVector statesWithProbabilityGreater0; + if (this->minimumOperatorStack.top()) { + statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound()); + } else { + statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound()); + } + + // Check if we already know the result (i.e. probability 0) for all initial states and + // don't compute anything in this case. + if (this->getInitialStates().isDisjointFrom(statesWithProbabilityGreater0)) { + LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step." + << " No exact probabilities were computed."); + // Set the values for all maybe-states to 0.5 to indicate that their probability values are not 0 (and + // not necessarily 1). + storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, Type(0.5)); + } else { + // In this case we have have to compute the probabilities. + + // We can eliminate the rows and columns from the original transition probability matrix that have probability 0. + storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0, this->getModel().getNondeterministicChoiceIndices()); + + // Get the "new" nondeterministic choice indices for the submatrix. + std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(statesWithProbabilityGreater0); + + // Compute the new set of target states in the reduced system. + storm::storage::BitVector rightStatesInReducedSystem = statesWithProbabilityGreater0 % *rightStates; + + // Make all rows absorbing that satisfy the second sub-formula. + submatrix.makeRowsAbsorbing(rightStatesInReducedSystem, subNondeterministicChoiceIndices); + + // Create the vector with which to multiply. + std::vector subresult(statesWithProbabilityGreater0.getNumberOfSetBits()); + storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne()); + + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->minimumOperatorStack.top(), submatrix, subresult, subNondeterministicChoiceIndices, nullptr, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Set the values of the resulting vector accordingly. + storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, subresult); + storm::utility::vector::setVectorValues(*result, ~statesWithProbabilityGreater0, storm::utility::constGetZero()); + } + + // Delete intermediate results and return result. + delete leftStates; + delete rightStates; + return result; + } + + /*! + * Checks the given formula that is a next formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @returns The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkNext(const storm::property::prctl::Next& formula, bool qualitative) const { + // First, we need to compute the states that satisfy the sub-formula of the next-formula. + storm::storage::BitVector* nextStates = formula.getChild().check(*this); + + // Create the vector with which to multiply and initialize it correctly. + std::vector* result = new std::vector(this->getModel().getNumberOfStates()); + storm::utility::vector::setVectorValues(*result, *nextStates, storm::utility::constGetOne()); + + // Delete obsolete sub-result. + delete nextStates; + + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->minimumOperatorStack.top(), this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Return result. + return result; + } + + /*! + * Checks the given formula that is a bounded-eventually formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @returns The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkBoundedEventually(const storm::property::prctl::BoundedEventually& formula, bool qualitative) const { + // Create equivalent temporary bounded until formula and check it. + storm::property::prctl::BoundedUntil temporaryBoundedUntilFormula(new storm::property::prctl::Ap("true"), formula.getChild().clone(), formula.getBound()); + return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative); + } + + /*! + * Checks the given formula that is an eventually formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @returns The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkEventually(const storm::property::prctl::Eventually& formula, bool qualitative) const { + // Create equivalent temporary until formula and check it. + storm::property::prctl::Until temporaryUntilFormula(new storm::property::prctl::Ap("true"), formula.getChild().clone()); + return this->checkUntil(temporaryUntilFormula, qualitative); + } + + /*! + * Checks the given formula that is a globally formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @returns The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkGlobally(const storm::property::prctl::Globally& formula, bool qualitative) const { + // Create "equivalent" temporary eventually formula and check it. + storm::property::prctl::Eventually temporaryEventuallyFormula(new storm::property::prctl::Not(formula.getChild().clone())); + std::vector* result = this->checkEventually(temporaryEventuallyFormula, qualitative); + + // Now subtract the resulting vector from the constant one vector to obtain final result. + storm::utility::vector::subtractFromConstantOneVector(*result); + return result; + } + + /*! + * Check the given formula that is an until formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @param scheduler If qualitative is false and this vector is non-null and has as many elements as + * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability + * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. + * @return The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkUntil(const storm::property::prctl::Until& formula, bool qualitative) const { + return this->checkUntil(this->minimumOperatorStack.top(), formula, qualitative, nullptr); + } + + /*! + * Check the given formula that is an until formula. + * + * @param minimize If set, the probability is minimized and maximized otherwise. + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bounds 0 and 1. + * @param scheduler If qualitative is false and this vector is non-null and has as many elements as + * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability + * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. + * @return The probabilities for the given formula to hold on every state of the model associated with this model + * checker. If the qualitative flag is set, exact probabilities might not be computed. + */ + virtual std::vector* checkUntil(bool minimize, const storm::property::prctl::Until& formula, bool qualitative, std::vector* scheduler) const { + // First, we need to compute the states that satisfy the sub-formulas of the until-formula. + storm::storage::BitVector* leftStates = formula.getLeft().check(*this); + storm::storage::BitVector* rightStates = formula.getRight().check(*this); + + // Then, we need to identify the states which have to be taken out of the matrix, i.e. + // all states that have probability 0 and 1 of satisfying the until-formula. + std::pair statesWithProbability01; + if (minimize) { + statesWithProbability01 = storm::utility::graph::performProb01Min(this->getModel(), *leftStates, *rightStates); + } else { + statesWithProbability01 = storm::utility::graph::performProb01Max(this->getModel(), *leftStates, *rightStates); + } + storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first); + storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second); + + // Delete sub-results that are obsolete now. + delete leftStates; + delete rightStates; + + storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1); + LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states."); + LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states."); + LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); + + // Create resulting vector. + std::vector* result = new std::vector(this->getModel().getNumberOfStates()); + + // Check whether we need to compute exact probabilities for some states. + if (this->getInitialStates().isDisjointFrom(maybeStates) || qualitative) { + if (qualitative) { + LOG4CPLUS_INFO(logger, "The formula was checked qualitatively. No exact probabilities were computed."); + } else { + LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step." + << " No exact probabilities were computed."); + } + // Set the values for all maybe-states to 0.5 to indicate that their probability values + // are neither 0 nor 1. + storm::utility::vector::setVectorValues(*result, maybeStates, Type(0.5)); + } else { + // In this case we have have to compute the probabilities. + + // First, we can eliminate the rows and columns from the original transition probability matrix for states + // whose probabilities are already known. + storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices()); + + // Get the "new" nondeterministic choice indices for the submatrix. + std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates); + + // Prepare the right-hand side of the equation system. For entry i this corresponds to + // the accumulated probability of going from state i to some 'yes' state. + 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(submatrix, subNondeterministicChoiceIndices, b); + + // Solve the corresponding system of equations. + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Set values of resulting vector according to result. + storm::utility::vector::setVectorValues(*result, maybeStates, x); + } + + // Set values of resulting vector that are known exactly. + storm::utility::vector::setVectorValues(*result, statesWithProbability0, storm::utility::constGetZero()); + storm::utility::vector::setVectorValues(*result, statesWithProbability1, storm::utility::constGetOne()); + + // If we were required to generate a scheduler, do so now. + if (scheduler != nullptr) { + this->computeTakenChoices(this->minimumOperatorStack.top(), *result, *scheduler, this->getModel().getNondeterministicChoiceIndices()); + } + + return result; + } + + /*! + * Checks the given formula that is an instantaneous reward formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bound 0. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bound 0. + * @returns The reward values for the given formula for every state of the model associated with this model + * checker. If the qualitative flag is set, exact values might not be computed. + */ + virtual std::vector* checkInstantaneousReward(const storm::property::prctl::InstantaneousReward& formula, bool qualitative) const { + // Only compute the result if the model has a state-based reward model. + if (!this->getModel().hasStateRewards()) { + LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula."); + throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula."; + } + + // Initialize result to state rewards of the model. + std::vector* result = new std::vector(this->getModel().getStateRewardVector()); + + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->minimumOperatorStack.top(), this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), nullptr, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Return result. + return result; + } + + /*! + * Check the given formula that is a cumulative reward formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bound 0. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bound 0. + * @returns The reward values for the given formula for every state of the model associated with this model + * checker. If the qualitative flag is set, exact values might not be computed. + */ + virtual std::vector* checkCumulativeReward(const storm::property::prctl::CumulativeReward& formula, bool qualitative) const { + // Only compute the result if the model has at least one reward model. + if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) { + LOG4CPLUS_ERROR(logger, "Missing reward model for formula."); + throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula."; + } + + // Compute the reward vector to add in each step based on the available reward models. + std::vector totalRewardVector; + if (this->getModel().hasTransitionRewards()) { + totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix()); + if (this->getModel().hasStateRewards()) { + storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector()); + } + } else { + totalRewardVector = std::vector(this->getModel().getStateRewardVector()); + } + + // Initialize result to either the state rewards of the model or the null vector. + std::vector* result = nullptr; + if (this->getModel().hasStateRewards()) { + result = new std::vector(this->getModel().getStateRewardVector()); + } else { + result = new std::vector(this->getModel().getNumberOfStates()); + } + + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->performMatrixVectorMultiplication(this->minimumOperatorStack.top(), this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), &totalRewardVector, formula.getBound()); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Delete temporary variables and return result. + return result; + } + + /*! + * Checks the given formula that is a reachability reward formula. + * + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bound 0. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bound 0. + * @return The reward values for the given formula for every state of the model associated with this model + * checker. If the qualitative flag is set, exact values might not be computed. + */ + virtual std::vector* checkReachabilityReward(const storm::property::prctl::ReachabilityReward& formula, bool qualitative) const { + return this->checkReachabilityReward(this->minimumOperatorStack.top(), formula, qualitative, nullptr); + } + + /*! + * Checks the given formula that is a reachability reward formula. + * + * @param minimize If set, the reward is to be minimized and maximized otherwise. + * @param formula The formula to check. + * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the + * results are only compared against the bound 0. If set to true, this will most likely results that are only + * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the + * bound 0. + * @param scheduler If qualitative is false and this vector is non-null and has as many elements as + * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability + * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. + * @return The reward values for the given formula for every state of the model associated with this model + * checker. If the qualitative flag is set, exact values might not be computed. + */ + virtual std::vector* checkReachabilityReward(bool minimize, const storm::property::prctl::ReachabilityReward& formula, bool qualitative, std::vector* scheduler) const { + // Only compute the result if the model has at least one reward model. + if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) { + LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula"); + throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula."; + } + + // Determine the states for which the target predicate holds. + storm::storage::BitVector* targetStates = formula.getChild().check(*this); + + // Determine which states have a reward of infinity by definition. + storm::storage::BitVector infinityStates; + storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true); + if (minimize) { + infinityStates = std::move(storm::utility::graph::performProb1A(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, *targetStates)); + } else { + infinityStates = std::move(storm::utility::graph::performProb1E(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, *targetStates)); + } + infinityStates.complement(); + storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates; + LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states."); + LOG4CPLUS_INFO(logger, "Found " << targetStates->getNumberOfSetBits() << " 'target' states."); + LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); + + // Create resulting vector. + std::vector* result = new std::vector(this->getModel().getNumberOfStates()); + + // Check whether we need to compute exact rewards for some states. + if (this->getInitialStates().isDisjointFrom(maybeStates)) { + LOG4CPLUS_INFO(logger, "The rewards for the initial states were determined in a preprocessing step." + << " No exact rewards were computed."); + // Set the values for all maybe-states to 1 to indicate that their reward values + // are neither 0 nor infinity. + storm::utility::vector::setVectorValues(*result, maybeStates, storm::utility::constGetOne()); + } else { + // In this case we have to compute the reward values for the remaining states. + + // We can eliminate the rows and columns from the original transition probability matrix for states + // whose reward values are already known. + storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices()); + + // Get the "new" nondeterministic choice indices for the submatrix. + std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates); + + // Prepare the right-hand side of the equation system. For entry i this corresponds to + // the accumulated probability of going from state i to some 'yes' state. + std::vector b(submatrix.getRowCount()); + + if (this->getModel().hasTransitionRewards()) { + // If a transition-based reward model is available, we initialize the right-hand + // side to the vector resulting from summing the rows of the pointwise product + // of the transition probability matrix and the transition reward matrix. + std::vector pointwiseProductRowSumVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix()); + storm::utility::vector::selectVectorValues(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), pointwiseProductRowSumVector); + + if (this->getModel().hasStateRewards()) { + // If a state-based reward model is also available, we need to add this vector + // as well. As the state reward vector contains entries not just for the states + // that we still consider (i.e. maybeStates), we need to extract these values + // first. + std::vector subStateRewards(b.size()); + storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector()); + storm::utility::vector::addVectorsInPlace(b, subStateRewards); + } + } else { + // If only a state-based reward model is available, we take this vector as the + // right-hand side. As the state reward vector contains entries not just for the + // states that we still consider (i.e. maybeStates), we need to extract these values + // first. + storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector()); + } + + // Create vector for results for maybe states. + std::vector x = this->getInitialValueIterationValues(submatrix, subNondeterministicChoiceIndices, b); + + // Solve the corresponding system of equations. + if (linearEquationSolver != nullptr) { + this->linearEquationSolver->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices); + } else { + throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available."; + } + + // Set values of resulting vector according to result. + storm::utility::vector::setVectorValues(*result, maybeStates, x); + } + + // Set values of resulting vector that are known exactly. + storm::utility::vector::setVectorValues(*result, *targetStates, storm::utility::constGetZero()); + storm::utility::vector::setVectorValues(*result, infinityStates, storm::utility::constGetInfinity()); + + // If we were required to generate a scheduler, do so now. + if (scheduler != nullptr) { + this->computeTakenChoices(this->minimumOperatorStack.top(), *result, *scheduler, this->getModel().getNondeterministicChoiceIndices()); + } + + // Delete temporary storages and return result. + delete targetStates; + return result; + } + + protected: + /*! + * Computes the vector of choices that need to be made to minimize/maximize the model checking result for each state. + * + * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. + * @param nondeterministicResult The model checking result for nondeterministic choices of all states. + * @param takenChoices The output vector that is to store the taken choices. + * @param nondeterministicChoiceIndices The assignment of states to their nondeterministic choices in the matrix. + */ + void computeTakenChoices(bool minimize, std::vector const& nondeterministicResult, std::vector& takenChoices, std::vector const& nondeterministicChoiceIndices) const { + std::vector temporaryResult(nondeterministicChoiceIndices.size() - 1); + if (minimize) { + storm::utility::vector::reduceVectorMin(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices); + } else { + storm::utility::vector::reduceVectorMax(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices); + } + } + + /*! + * A stack used for storing whether we are currently computing min or max probabilities or rewards, respectively. + * The topmost element is true if and only if we are currently computing minimum probabilities or rewards. + */ + mutable std::stack minimumOperatorStack; + + private: + + /*! + * Computes the nondeterministic choice indices vector resulting from reducing the full system to the states given + * by the parameter constraint. + * + * @param constraint A bit vector specifying which states are kept. + * @returns A vector of the nondeterministic choice indices of the subsystem induced by the given constraint. + */ + std::vector computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector const& constraint) const { + // First, get a reference to the full nondeterministic choice indices. + std::vector const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); + + // Reserve the known amount of slots for the resulting vector. + std::vector subNondeterministicChoiceIndices(constraint.getNumberOfSetBits() + 1); + uint_fast64_t currentRowCount = 0; + uint_fast64_t currentIndexCount = 1; + + // Set the first element as this will clearly begin at offset 0. + subNondeterministicChoiceIndices[0] = 0; + + // Loop over all states that need to be kept and copy the relative indices of the nondeterministic choices over + // to the resulting vector. + for (auto index : constraint) { + subNondeterministicChoiceIndices[currentIndexCount] = currentRowCount + nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index]; + currentRowCount += nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index]; + ++currentIndexCount; + } + + // Put a sentinel element at the end. + subNondeterministicChoiceIndices[constraint.getNumberOfSetBits()] = currentRowCount; + + return subNondeterministicChoiceIndices; + } + + /*! + * Retrieves the values to be used as the initial values for the value iteration techniques. + * + * @param submatrix The matrix that will be used for value iteration later. + */ + std::vector getInitialValueIterationValues(storm::storage::SparseMatrix const& submatrix, std::vector const& subNondeterministicChoiceIndices, std::vector const& rightHandSide) const { + // std::vector scheduler(submatrix.getColumnCount()); + // std::vector result(scheduler.size(), Type(0.5)); + // std::vector b(scheduler.size()); + // storm::utility::vector::selectVectorValues(b, scheduler, subNondeterministicChoiceIndices, rightHandSide); + // storm::storage::SparseMatrix A(submatrix.getSubmatrix(scheduler, subNondeterministicChoiceIndices)); + // A.convertToEquationSystem(); + // this->solveEquationSystem(A, result, b); + std::vector result(submatrix.getColumnCount()); + return result; + } + + // An object that is used for solving linear equations and performing matrix-vector multiplication. + std::unique_ptr> linearEquationSolver; + + }; - // Create the vector with which to multiply. - std::vector subresult(statesWithProbabilityGreater0.getNumberOfSetBits()); - storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne()); - - this->performMatrixVectorMultiplication(submatrix, subresult, subNondeterministicChoiceIndices, nullptr, formula.getBound()); - - // Set the values of the resulting vector accordingly. - storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, subresult); - storm::utility::vector::setVectorValues(*result, ~statesWithProbabilityGreater0, storm::utility::constGetZero()); - } - - // Delete intermediate results and return result. - delete leftStates; - delete rightStates; - return result; - } - - /*! - * Checks the given formula that is a next formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @returns The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkNext(const storm::property::prctl::Next& formula, bool qualitative) const { - // First, we need to compute the states that satisfy the sub-formula of the next-formula. - storm::storage::BitVector* nextStates = formula.getChild().check(*this); - - // Create the vector with which to multiply and initialize it correctly. - std::vector* result = new std::vector(this->getModel().getNumberOfStates()); - storm::utility::vector::setVectorValues(*result, *nextStates, storm::utility::constGetOne()); - - // Delete obsolete sub-result. - delete nextStates; - - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices()); - - // Return result. - return result; - } - - /*! - * Checks the given formula that is a bounded-eventually formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @returns The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkBoundedEventually(const storm::property::prctl::BoundedEventually& formula, bool qualitative) const { - // Create equivalent temporary bounded until formula and check it. - storm::property::prctl::BoundedUntil temporaryBoundedUntilFormula(new storm::property::prctl::Ap("true"), formula.getChild().clone(), formula.getBound()); - return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative); - } - - /*! - * Checks the given formula that is an eventually formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @returns The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkEventually(const storm::property::prctl::Eventually& formula, bool qualitative) const { - // Create equivalent temporary until formula and check it. - storm::property::prctl::Until temporaryUntilFormula(new storm::property::prctl::Ap("true"), formula.getChild().clone()); - return this->checkUntil(temporaryUntilFormula, qualitative); - } - - /*! - * Checks the given formula that is a globally formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @returns The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkGlobally(const storm::property::prctl::Globally& formula, bool qualitative) const { - // Create "equivalent" temporary eventually formula and check it. - storm::property::prctl::Eventually temporaryEventuallyFormula(new storm::property::prctl::Not(formula.getChild().clone())); - std::vector* result = this->checkEventually(temporaryEventuallyFormula, qualitative); - - // Now subtract the resulting vector from the constant one vector to obtain final result. - storm::utility::vector::subtractFromConstantOneVector(*result); - return result; - } - - /*! - * Check the given formula that is an until formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @param scheduler If qualitative is false and this vector is non-null and has as many elements as - * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability - * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. - * @return The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkUntil(const storm::property::prctl::Until& formula, bool qualitative) const { - return this->checkUntil(this->minimumOperatorStack.top(), formula, qualitative, nullptr); - } - - /*! - * Check the given formula that is an until formula. - * - * @param minimize If set, the probability is minimized and maximized otherwise. - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bounds 0 and 1. - * @param scheduler If qualitative is false and this vector is non-null and has as many elements as - * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability - * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. - * @return The probabilities for the given formula to hold on every state of the model associated with this model - * checker. If the qualitative flag is set, exact probabilities might not be computed. - */ - virtual std::vector* checkUntil(bool minimize, const storm::property::prctl::Until& formula, bool qualitative, std::vector* scheduler) const { - // First, we need to compute the states that satisfy the sub-formulas of the until-formula. - storm::storage::BitVector* leftStates = formula.getLeft().check(*this); - storm::storage::BitVector* rightStates = formula.getRight().check(*this); - - // Then, we need to identify the states which have to be taken out of the matrix, i.e. - // all states that have probability 0 and 1 of satisfying the until-formula. - std::pair statesWithProbability01; - if (minimize) { - statesWithProbability01 = storm::utility::graph::performProb01Min(this->getModel(), *leftStates, *rightStates); - } else { - statesWithProbability01 = storm::utility::graph::performProb01Max(this->getModel(), *leftStates, *rightStates); - } - storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first); - storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second); - - // Delete sub-results that are obsolete now. - delete leftStates; - delete rightStates; - - storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1); - LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states."); - LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states."); - LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); - - // Create resulting vector. - std::vector* result = new std::vector(this->getModel().getNumberOfStates()); - - // Check whether we need to compute exact probabilities for some states. - if (this->getInitialStates().isDisjointFrom(maybeStates) || qualitative) { - if (qualitative) { - LOG4CPLUS_INFO(logger, "The formula was checked qualitatively. No exact probabilities were computed."); - } else { - LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step." - << " No exact probabilities were computed."); - } - // Set the values for all maybe-states to 0.5 to indicate that their probability values - // are neither 0 nor 1. - storm::utility::vector::setVectorValues(*result, maybeStates, Type(0.5)); - } else { - // In this case we have have to compute the probabilities. - - // First, we can eliminate the rows and columns from the original transition probability matrix for states - // whose probabilities are already known. - storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices()); - - // Get the "new" nondeterministic choice indices for the submatrix. - std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates); - - // Prepare the right-hand side of the equation system. For entry i this corresponds to - // the accumulated probability of going from state i to some 'yes' state. - 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(submatrix, formula); - - // Solve the corresponding system of equations. - this->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices); - - // Set values of resulting vector according to result. - storm::utility::vector::setVectorValues(*result, maybeStates, x); - } - - // Set values of resulting vector that are known exactly. - storm::utility::vector::setVectorValues(*result, statesWithProbability0, storm::utility::constGetZero()); - storm::utility::vector::setVectorValues(*result, statesWithProbability1, storm::utility::constGetOne()); - - return result; - } - - /*! - * Checks the given formula that is an instantaneous reward formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bound 0. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bound 0. - * @returns The reward values for the given formula for every state of the model associated with this model - * checker. If the qualitative flag is set, exact values might not be computed. - */ - virtual std::vector* checkInstantaneousReward(const storm::property::prctl::InstantaneousReward& formula, bool qualitative) const { - // Only compute the result if the model has a state-based reward model. - if (!this->getModel().hasStateRewards()) { - LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula."); - throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula."; - } - - // Initialize result to state rewards of the model. - std::vector* result = new std::vector(this->getModel().getStateRewardVector()); - - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), nullptr, formula.getBound()); - - // Return result. - return result; - } - - /*! - * Check the given formula that is a cumulative reward formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bound 0. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bound 0. - * @returns The reward values for the given formula for every state of the model associated with this model - * checker. If the qualitative flag is set, exact values might not be computed. - */ - virtual std::vector* checkCumulativeReward(const storm::property::prctl::CumulativeReward& formula, bool qualitative) const { - // Only compute the result if the model has at least one reward model. - if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) { - LOG4CPLUS_ERROR(logger, "Missing reward model for formula."); - throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula."; - } - - // Compute the reward vector to add in each step based on the available reward models. - std::vector totalRewardVector; - if (this->getModel().hasTransitionRewards()) { - totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix()); - if (this->getModel().hasStateRewards()) { - storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector()); - } - } else { - totalRewardVector = std::vector(this->getModel().getStateRewardVector()); - } - - // Initialize result to either the state rewards of the model or the null vector. - std::vector* result = nullptr; - if (this->getModel().hasStateRewards()) { - result = new std::vector(this->getModel().getStateRewardVector()); - } else { - result = new std::vector(this->getModel().getNumberOfStates()); - } - - this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), &totalRewardVector, formula.getBound()); - - // Delete temporary variables and return result. - return result; - } - - /*! - * Checks the given formula that is a reachability reward formula. - * - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bound 0. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bound 0. - * @return The reward values for the given formula for every state of the model associated with this model - * checker. If the qualitative flag is set, exact values might not be computed. - */ - virtual std::vector* checkReachabilityReward(const storm::property::prctl::ReachabilityReward& formula, bool qualitative) const { - return this->checkReachabilityReward(this->minimumOperatorStack.top(), formula, qualitative, nullptr); - } - - /*! - * Checks the given formula that is a reachability reward formula. - * - * @param minimize If set, the reward is to be minimized and maximized otherwise. - * @param formula The formula to check. - * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the - * results are only compared against the bound 0. If set to true, this will most likely results that are only - * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the - * bound 0. - * @param scheduler If qualitative is false and this vector is non-null and has as many elements as - * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability - * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state. - * @return The reward values for the given formula for every state of the model associated with this model - * checker. If the qualitative flag is set, exact values might not be computed. - */ - virtual std::vector* checkReachabilityReward(bool minimize, const storm::property::prctl::ReachabilityReward& formula, bool qualitative, std::vector* scheduler) const { - // Only compute the result if the model has at least one reward model. - if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) { - LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula"); - throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula."; - } - - // Determine the states for which the target predicate holds. - storm::storage::BitVector* targetStates = formula.getChild().check(*this); - - // Determine which states have a reward of infinity by definition. - storm::storage::BitVector infinityStates; - storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true); - if (minimize) { - infinityStates = std::move(storm::utility::graph::performProb1A(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, *targetStates)); - } else { - infinityStates = std::move(storm::utility::graph::performProb1E(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, *targetStates)); - } - infinityStates.complement(); - storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates; - LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states."); - LOG4CPLUS_INFO(logger, "Found " << targetStates->getNumberOfSetBits() << " 'target' states."); - LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); - - // Create resulting vector. - std::vector* result = new std::vector(this->getModel().getNumberOfStates()); - - // Check whether we need to compute exact rewards for some states. - if (this->getInitialStates().isDisjointFrom(maybeStates)) { - LOG4CPLUS_INFO(logger, "The rewards for the initial states were determined in a preprocessing step." - << " No exact rewards were computed."); - // Set the values for all maybe-states to 1 to indicate that their reward values - // are neither 0 nor infinity. - storm::utility::vector::setVectorValues(*result, maybeStates, storm::utility::constGetOne()); - } else { - // In this case we have to compute the reward values for the remaining states. - - // We can eliminate the rows and columns from the original transition probability matrix for states - // whose reward values are already known. - storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices()); - - // Get the "new" nondeterministic choice indices for the submatrix. - std::vector subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates); - - // Prepare the right-hand side of the equation system. For entry i this corresponds to - // the accumulated probability of going from state i to some 'yes' state. - std::vector b(submatrix.getRowCount()); - - if (this->getModel().hasTransitionRewards()) { - // If a transition-based reward model is available, we initialize the right-hand - // side to the vector resulting from summing the rows of the pointwise product - // of the transition probability matrix and the transition reward matrix. - std::vector pointwiseProductRowSumVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix()); - storm::utility::vector::selectVectorValues(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), pointwiseProductRowSumVector); - - if (this->getModel().hasStateRewards()) { - // If a state-based reward model is also available, we need to add this vector - // as well. As the state reward vector contains entries not just for the states - // that we still consider (i.e. maybeStates), we need to extract these values - // first. - std::vector subStateRewards(b.size()); - storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector()); - storm::utility::vector::addVectorsInPlace(b, subStateRewards); - } - } else { - // If only a state-based reward model is available, we take this vector as the - // right-hand side. As the state reward vector contains entries not just for the - // states that we still consider (i.e. maybeStates), we need to extract these values - // first. - storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector()); - } - - // Create vector for results for maybe states. - std::vector x = this->getInitialValueIterationValues(submatrix, formula); - - // Solve the corresponding system of equations. - this->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices, scheduler); - - // Set values of resulting vector according to result. - storm::utility::vector::setVectorValues(*result, maybeStates, x); - } - - // Set values of resulting vector that are known exactly. - storm::utility::vector::setVectorValues(*result, *targetStates, storm::utility::constGetZero()); - storm::utility::vector::setVectorValues(*result, infinityStates, storm::utility::constGetInfinity()); - - // Delete temporary storages and return result. - delete targetStates; - return result; - } - -protected: - /*! - * Computes the vector of choices that need to be made to minimize/maximize the model checking result for each state. - * - * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. - * @param nondeterministicResult The model checking result for nondeterministic choices of all states. - * @param takenChoices The output vector that is to store the taken choices. - * @param nondeterministicChoiceIndices The assignment of states to their nondeterministic choices in the matrix. - */ - void computeTakenChoices(bool minimize, std::vector const& nondeterministicResult, std::vector& takenChoices, std::vector const& nondeterministicChoiceIndices) const { - std::vector temporaryResult(nondeterministicChoiceIndices.size() - 1); - if (minimize) { - storm::utility::vector::reduceVectorMin(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices); - } else { - storm::utility::vector::reduceVectorMax(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices); - } - } - - /*! - * A stack used for storing whether we are currently computing min or max probabilities or rewards, respectively. - * The topmost element is true if and only if we are currently computing minimum probabilities or rewards. - */ - mutable std::stack minimumOperatorStack; - -private: - /*! - * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b - * until x[n], where x[0] = x. - * - * @param A The matrix that is to be multiplied against the vector. - * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, - * i.e. after the method returns, this vector will contain the computed values. - * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. - * @param b If not null, this vector is being added to the result after each matrix-vector multiplication. - * @param n Specifies the number of iterations the matrix-vector multiplication is performed. - * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. - */ - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& nondeterministicChoiceIndices, std::vector* b = nullptr, uint_fast64_t n = 1) const { - // Create vector for result of multiplication, which is reduced to the result vector after - // each multiplication. - std::vector multiplyResult(A.getRowCount()); - - // Now perform matrix-vector multiplication as long as we meet the bound of the formula. - for (uint_fast64_t i = 0; i < n; ++i) { - A.multiplyWithVector(x, multiplyResult); - - // Add b if it is non-null. - if (b != nullptr) { - storm::utility::vector::addVectorsInPlace(multiplyResult, *b); - } - - // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost - // element of the min/max operator stack. - if (this->minimumOperatorStack.top()) { - storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); - } else { - storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); - } - } - } - - /*! - * Solves the equation system A*x = b given by the parameters. - * - * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. - * @param A The matrix specifying the coefficients of the linear equations. - * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but - * may be ignored. - * @param b The right-hand side of the equation system. - * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. - * @param takenChoices If not null, this vector will be filled with the nondeterministic choices taken by the states - * to achieve the solution of the equation system. This assumes that the given vector has at least as many elements - * as there are states in the MDP. - * @returns The solution vector x of the system of linear equations as the content of the parameter x. - */ - virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector const& nondeterministicChoiceIndices, std::vector* takenChoices = nullptr) const { - LOG4CPLUS_INFO(logger, "Starting iterative solver."); - - // Get the settings object to customize solving. - storm::settings::Settings* s = storm::settings::instance(); - - // Get relevant user-defined settings for solving the equations. - double precision = s->get("precision"); - unsigned maxIterations = s->get("maxiter"); - bool relative = s->get("relative"); - - // Set up the environment for the power method. - std::vector multiplyResult(A.getRowCount()); - std::vector* currentX = &x; - std::vector* newX = new std::vector(x.size()); - std::vector* swap = nullptr; - uint_fast64_t iterations = 0; - bool converged = false; - - // Proceed with the iterations as long as the method did not converge or reach the - // user-specified maximum number of iterations. - while (!converged && iterations < maxIterations) { - // Compute x' = A*x + b. - A.multiplyWithVector(*currentX, multiplyResult); - storm::utility::vector::addVectorsInPlace(multiplyResult, b); - - // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost - // element of the min/max operator stack. - if (minimize) { - storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); - } else { - storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); - } - - // Determine whether the method converged. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); - - // Update environment variables. - swap = currentX; - currentX = newX; - newX = swap; - ++iterations; - } - - // If we were requested to record the taken choices, we have to construct the vector now. - if (takenChoices != nullptr) { - this->computeTakenChoices(minimize, multiplyResult, *takenChoices, nondeterministicChoiceIndices); - } - - // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result - // is currently stored in currentX, but x is the output vector. - if (iterations % 2 == 1) { - std::swap(x, *currentX); - delete currentX; - } else { - delete newX; - } - - // Check if the solver converged and issue a warning otherwise. - if (converged) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - } - - /*! - * Computes the nondeterministic choice indices vector resulting from reducing the full system to the states given - * by the parameter constraint. - * - * @param constraint A bit vector specifying which states are kept. - * @returns A vector of the nondeterministic choice indices of the subsystem induced by the given constraint. - */ - std::vector computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector const& constraint) const { - // First, get a reference to the full nondeterministic choice indices. - std::vector const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); - - // Reserve the known amount of slots for the resulting vector. - std::vector subNondeterministicChoiceIndices(constraint.getNumberOfSetBits() + 1); - uint_fast64_t currentRowCount = 0; - uint_fast64_t currentIndexCount = 1; - - // Set the first element as this will clearly begin at offset 0. - subNondeterministicChoiceIndices[0] = 0; - - // Loop over all states that need to be kept and copy the relative indices of the nondeterministic choices over - // to the resulting vector. - for (auto index : constraint) { - subNondeterministicChoiceIndices[currentIndexCount] = currentRowCount + nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index]; - currentRowCount += nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index]; - ++currentIndexCount; - } - - // Put a sentinel element at the end. - subNondeterministicChoiceIndices[constraint.getNumberOfSetBits()] = currentRowCount; - - return subNondeterministicChoiceIndices; - } - - /*! - * Retrieves the values to be used as the initial values for the value iteration techniques. - * - * @param submatrix The matrix that will be used for value iteration later. - */ - std::vector getInitialValueIterationValues(storm::storage::SparseMatrix const& submatrix, storm::property::prctl::AbstractPathFormula const& formula) const { - // std::vector scheduler(submatrix.getColumnCount()); - std::vector scheduler(this->getModel().getNumberOfStates()); - - storm::models::Dtmc dtmc(this->getModel().getTransitionMatrix().getSubmatrix(scheduler, this->getModel().getNondeterministicChoiceIndices()), - storm::models::AtomicPropositionsLabeling(this->getModel().getStateLabeling()), - this->getModel().hasStateRewards() ? boost::optional>(this->getModel().getStateRewardVector()) : boost::optional>(), - this->getModel().hasTransitionRewards() ? boost::optional>(this->getModel().getTransitionRewardMatrix().getSubmatrix(scheduler, this->getModel().getNondeterministicChoiceIndices(), false)) : boost::optional>()); - - storm::modelchecker::prctl::GmmxxDtmcPrctlModelChecker modelChecker(dtmc); - - std::vector* modelCheckingResult = formula.check(modelChecker, false); - std::vector result(std::move(*modelCheckingResult)); - delete modelCheckingResult; - - return result; - } - -}; - -} // namespace prctl -} // namespace modelchecker + } // namespace prctl + } // namespace modelchecker } // namespace storm #endif /* STORM_MODELCHECKER_PRCTL_SPARSEMDPPRCTLMODELCHECKER_H_ */ diff --git a/src/models/AtomicPropositionsLabeling.h b/src/models/AtomicPropositionsLabeling.h index e0e423be7..e566dfe7d 100644 --- a/src/models/AtomicPropositionsLabeling.h +++ b/src/models/AtomicPropositionsLabeling.h @@ -58,6 +58,23 @@ public: singleLabelings(atomicPropositionsLabeling.singleLabelings) { // Intentionally left empty. } + + /*! + * Copy constructor that performs a deep copy of the given atomic proposition labeling while restricting to the states + * given by the bit vector. + * + * @param atomicPropositionsLabeling The atomic propositions labeling to copy. + * @param substates A subset of the states that is to be copied. + */ + AtomicPropositionsLabeling(AtomicPropositionsLabeling const& atomicPropositionsLabeling, storm::storage::BitVector const& substates) + : stateCount(substates.getNumberOfSetBits()), apCountMax(atomicPropositionsLabeling.apCountMax), apsCurrent(atomicPropositionsLabeling.apsCurrent), + nameToLabelingMap(atomicPropositionsLabeling.nameToLabelingMap) { + // Now we need to copy the fraction of the single labelings given by the substates. + singleLabelings.reserve(apCountMax); + for (auto const& labeling : atomicPropositionsLabeling.singleLabelings) { + singleLabelings.emplace_back(labeling, substates); + } + } /*! * Move Constructor that performs a move copy on the given atomic proposition labeling. diff --git a/src/solver/AbstractLinearEquationSolver.h b/src/solver/AbstractLinearEquationSolver.h new file mode 100644 index 000000000..25905e29e --- /dev/null +++ b/src/solver/AbstractLinearEquationSolver.h @@ -0,0 +1,26 @@ +#ifndef STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ + +#include "src/storage/SparseMatrix.h" + +#include + +namespace storm { + namespace solver { + + template + class AbstractLinearEquationSolver { + public: + + virtual AbstractLinearEquationSolver* clone() const = 0; + + virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const = 0; + + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1) const = 0; + + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ */ diff --git a/src/solver/AbstractNondeterministicLinearEquationSolver.h b/src/solver/AbstractNondeterministicLinearEquationSolver.h new file mode 100644 index 000000000..09ee98574 --- /dev/null +++ b/src/solver/AbstractNondeterministicLinearEquationSolver.h @@ -0,0 +1,175 @@ +#ifndef STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ + +#include "src/storage/SparseMatrix.h" + +#include + +namespace storm { + namespace solver { + + template + class AbstractNondeterministicLinearEquationSolver { + public: + + virtual AbstractNondeterministicLinearEquationSolver* clone() const { + return new AbstractNondeterministicLinearEquationSolver(); + } + + /*! + * Solves the equation system A*x = b given by the parameters. + * + * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. + * @param A The matrix specifying the coefficients of the linear equations. + * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but + * may be ignored. + * @param b The right-hand side of the equation system. + * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. + * @param takenChoices If not null, this vector will be filled with the nondeterministic choices taken by the states + * to achieve the solution of the equation system. This assumes that the given vector has at least as many elements + * as there are states in the MDP. + * @returns The solution vector x of the system of linear equations as the content of the parameter x. + */ + virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector const& nondeterministicChoiceIndices) const { + LOG4CPLUS_INFO(logger, "Starting iterative solver."); + + // Get the settings object to customize solving. + storm::settings::Settings* s = storm::settings::instance(); + + // Get relevant user-defined settings for solving the equations. + double precision = s->get("precision"); + unsigned maxIterations = s->get("maxiter"); + bool relative = s->get("relative"); + + // Set up the environment for the power method. + std::vector multiplyResult(A.getRowCount()); + std::vector* currentX = &x; + std::vector* newX = new std::vector(x.size()); + std::vector* swap = nullptr; + uint_fast64_t iterations = 0; + bool converged = false; + + // Proceed with the iterations as long as the method did not converge or reach the + // user-specified maximum number of iterations. + while (!converged && iterations < maxIterations) { + // Compute x' = A*x + b. + A.multiplyWithVector(*currentX, multiplyResult); + storm::utility::vector::addVectorsInPlace(multiplyResult, b); + + // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost + // element of the min/max operator stack. + if (minimize) { + storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); + } else { + storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); + + // Update environment variables. + swap = currentX; + currentX = newX; + newX = swap; + ++iterations; + } + + // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result + // is currently stored in currentX, but x is the output vector. + if (iterations % 2 == 1) { + std::swap(x, *currentX); + delete currentX; + } else { + delete newX; + } + + // Check if the solver converged and issue a warning otherwise. + if (converged) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } + + /*! + * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b + * until x[n], where x[0] = x. + * + * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. + * @param A The matrix that is to be multiplied against the vector. + * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, + * i.e. after the method returns, this vector will contain the computed values. + * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. + * @param b If not null, this vector is being added to the result after each matrix-vector multiplication. + * @param n Specifies the number of iterations the matrix-vector multiplication is performed. + * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. + */ + virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& nondeterministicChoiceIndices, std::vector* b = nullptr, uint_fast64_t n = 1) const { + // Create vector for result of multiplication, which is reduced to the result vector after + // each multiplication. + std::vector multiplyResult(A.getRowCount()); + + // Now perform matrix-vector multiplication as long as we meet the bound of the formula. + for (uint_fast64_t i = 0; i < n; ++i) { + A.multiplyWithVector(x, multiplyResult); + + // Add b if it is non-null. + if (b != nullptr) { + storm::utility::vector::addVectorsInPlace(multiplyResult, *b); + } + + // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost + // element of the min/max operator stack. + if (minimize) { + storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); + } else { + storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); + } + } + } + + /*! + * Returns the name of this solver. + * + * @returns The name of this solver. + */ + static std::string getName() { + return "native"; + } + + /*! + * Registers all options associated with the gmm++ matrix library. + */ + static void putOptions(boost::program_options::options_description* desc) { + desc->add_options()("lemethod", boost::program_options::value()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}."); + desc->add_options()("maxiter", boost::program_options::value()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); + desc->add_options()("precision", boost::program_options::value()->default_value(1e-6), "Sets the precision for iterative equation solving."); + desc->add_options()("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}."); + desc->add_options()("relative", boost::program_options::value()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence."); + } + + /*! + * Validates whether the given lemethod matches one of the available ones. + * 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")) { + throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid."; + } + } + + /*! + * Validates whether the given preconditioner matches one of the available ones. + * Throws an exception of type InvalidSettings in case the selected preconditioner is illegal. + */ + static void validatePreconditioner(const std::string& preconditioner) { + if ((preconditioner != "ilu") && (preconditioner != "diagonal") && (preconditioner != "ildlt") && (preconditioner != "none")) { + throw exceptions::InvalidSettingsException() << "Argument " << preconditioner << " for option 'precond' is invalid."; + } + } + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ */ diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h new file mode 100644 index 000000000..0c1570722 --- /dev/null +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -0,0 +1,228 @@ +#ifndef STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ + +#include "AbstractLinearEquationSolver.h" +#include "src/adapters/GmmxxAdapter.h" +#include "src/utility/ConstTemplates.h" +#include "src/utility/Settings.h" + +#include "gmm/gmm_matrix.h" +#include "gmm/gmm_iter_solvers.h" + +#include + +namespace storm { + namespace solver { + + template + class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver { + public: + + virtual AbstractLinearEquationSolver* clone() const { + return new GmmxxLinearEquationSolver(); + } + + virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { + // Get the settings object to customize linear solving. + storm::settings::Settings* s = storm::settings::instance(); + + // Prepare an iteration object that determines the accuracy, maximum number of iterations + // and the like. + gmm::iteration iter(s->get("precision"), 0, s->get("maxiter")); + + // Print some information about the used preconditioner. + const std::string& precond = s->getString("precond"); + LOG4CPLUS_INFO(logger, "Starting iterative solver."); + if (s->getString("lemethod") == "jacobi") { + if (precond != "none") { + LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the Jacobi method. Dropping preconditioner."); + } + } else { + if (precond == "ilu") { + LOG4CPLUS_INFO(logger, "Using ILU preconditioner."); + } else if (precond == "diagonal") { + LOG4CPLUS_INFO(logger, "Using diagonal preconditioner."); + } else if (precond == "ildlt") { + LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner."); + } else if (precond == "none") { + LOG4CPLUS_INFO(logger, "Using no preconditioner."); + } + } + + // Now do the actual solving. + if (s->getString("lemethod") == "bicgstab") { + LOG4CPLUS_INFO(logger, "Using BiCGStab 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::bicgstab(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); + } else if (precond == "diagonal") { + gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); + } else if (precond == "ildlt") { + gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); + } else if (precond == "none") { + gmm::bicgstab(*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") == "qmr") { + LOG4CPLUS_INFO(logger, "Using QMR 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::qmr(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); + } else if (precond == "diagonal") { + gmm::qmr(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); + } else if (precond == "ildlt") { + gmm::qmr(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); + } else if (precond == "none") { + 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") == "jacobi") { + LOG4CPLUS_INFO(logger, "Using Jacobi method."); + uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b); + + // Check if the solver converged and issue a warning otherwise. + if (iterations < s->get("maxiter")) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } + } + + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b, uint_fast64_t n = 1) const { + // Transform the transition probability A to the gmm++ format to use its arithmetic. + gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + + // Set up some temporary variables so that we can just swap pointers instead of copying the result after each + // iteration. + std::vector* swap = nullptr; + std::vector* currentVector = &x; + std::vector* tmpVector = new std::vector(A.getRowCount()); + + // Now perform matrix-vector multiplication as long as we meet the bound. + for (uint_fast64_t i = 0; i < n; ++i) { + gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector); + swap = tmpVector; + tmpVector = currentVector; + currentVector = swap; + + // If requested, add an offset to the current result vector. + if (b != nullptr) { + gmm::add(*b, *currentVector); + } + } + + // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, because + // the output is supposed to be stored in x. + if (n % 2 == 1) { + std::swap(x, *currentVector); + delete currentVector; + } else { + delete tmpVector; + } + + delete gmmxxMatrix; + } + + private: + /*! + * Solves the linear equation system A*x = b given by the parameters using the Jacobi method. + * + * @param A The matrix specifying the coefficients of the linear equations. + * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but + * may be ignored. + * @param b The right-hand side of the equation system. + * @returns The solution vector x of the system of linear equations as the content of the parameter x. + * @returns The number of iterations needed until convergence. + */ + uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { + // Get the settings object to customize linear solving. + storm::settings::Settings* s = storm::settings::instance(); + + double precision = s->get("precision"); + uint_fast64_t maxIterations = s->get("maxiter"); + bool relative = s->get("relative"); + + // Get a Jacobi decomposition of the matrix A. + typename storm::storage::SparseMatrix::SparseJacobiDecomposition_t jacobiDecomposition = A.getJacobiDecomposition(); + + // Convert the (inverted) diagonal matrix to gmm++'s format. + gmm::csr_matrix* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.second)); + // Convert the LU matrix to gmm++'s format. + gmm::csr_matrix* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); + + LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver."); + + // x_(k + 1) = D^-1 * (b - R * x_k) + // In x we keep a copy of the result for swapping in the loop (e.g. less copy-back). + std::vector* xNext = new std::vector(x.size()); + const std::vector* xCopy = xNext; + std::vector* xCurrent = &x; + + // Target vector for precision calculation. + std::vector* residuum = new std::vector(x.size()); + + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < maxIterations) { + // R * x_k (xCurrent is x_k) -> xNext + gmm::mult(*gmmxxLU, *xCurrent, *xNext); + // b - R * x_k (xNext contains R * x_k) -> xNext + gmm::add(b, gmm::scaled(*xNext, -1.0), *xNext); + // D^-1 * (b - R * x_k) -> xNext + gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext); + + // Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the + // vector. + std::vector *const swap = xNext; + xNext = xCurrent; + xCurrent = swap; + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(*xCurrent, *xNext, precision, relative); + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // 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 parameter x. + if (xCurrent == xCopy) { + std::swap(x, *xCurrent); + } + + // As xCopy always points to the copy of x used for swapping, we can safely delete it. + delete xCopy; + + // Also delete the other dynamically allocated variables. + delete residuum; + delete gmmxxDiagonalInverted; + delete gmmxxLU; + + return iterationCount; + } + + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ */ diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.h b/src/solver/GmmxxNondeterministicLinearEquationSolver.h new file mode 100644 index 000000000..e2cf5afb4 --- /dev/null +++ b/src/solver/GmmxxNondeterministicLinearEquationSolver.h @@ -0,0 +1,122 @@ +#ifndef STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ + +#include "AbstractLinearEquationSolver.h" +#include "src/adapters/GmmxxAdapter.h" +#include "src/storage/SparseMatrix.h" + +#include "gmm/gmm_matrix.h" + +#include + +namespace storm { + namespace solver { + + template + class GmmxxNondeterministicLinearEquationSolver : public AbstractNondeterministicLinearEquationSolver { + public: + + virtual AbstractLinearEquationSolver* clone() const { + return new GmmxxNondeterministicLinearEquationSolver(); + } + + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& nondeterministicChoiceIndices, std::vector* b = nullptr, uint_fast64_t n = 1) const { + // Transform the transition probability matrix to the gmm++ format to use its arithmetic. + gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + + // Create vector for result of multiplication, which is reduced to the result vector after + // each multiplication. + std::vector multiplyResult(A.getRowCount()); + + // Now perform matrix-vector multiplication as long as we meet the bound of the formula. + for (uint_fast64_t i = 0; i < n; ++i) { + gmm::mult(*gmmxxMatrix, x, multiplyResult); + + if (b != nullptr) { + gmm::add(*b, multiplyResult); + } + + if (this->minimumOperatorStack.top()) { + storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); + } else { + storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); + } + } + + // Delete intermediate results and return result. + delete gmmxxMatrix; + } + + void solveEquationSystem(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector const& nondeterministicChoiceIndices, std::vector* takenChoices = nullptr) const override { + // Get the settings object to customize solving. + storm::settings::Settings* s = storm::settings::instance(); + + // Get relevant user-defined settings for solving the equations. + double precision = s->get("precision"); + unsigned maxIterations = s->get("maxiter"); + bool relative = s->get("relative"); + + // Transform the transition probability matrix to the gmm++ format to use its arithmetic. + gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + + // Set up the environment for the power method. + std::vector multiplyResult(A.getRowCount()); + std::vector* currentX = &x; + std::vector* newX = new std::vector(x.size()); + std::vector* swap = nullptr; + uint_fast64_t iterations = 0; + bool converged = false; + + // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number + // of iterations. + while (!converged && iterations < maxIterations) { + // Compute x' = A*x + b. + gmm::mult(*gmmxxMatrix, *currentX, multiplyResult); + gmm::add(b, multiplyResult); + + // Reduce the vector x' by applying min/max for all non-deterministic choices. + if (minimize) { + storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); + } else { + storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); + + // Update environment variables. + swap = currentX; + currentX = newX; + newX = swap; + ++iterations; + } + + // If we were requested to record the taken choices, we have to construct the vector now. + if (takenChoices != nullptr) { + this->computeTakenChoices(minimize, multiplyResult, *takenChoices, nondeterministicChoiceIndices); + } + + // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result + // is currently stored in currentX, but x is the output vector. + if (iterations % 2 == 1) { + std::swap(x, *currentX); + delete currentX; + } else { + delete newX; + } + delete gmmxxMatrix; + + // Check if the solver converged and issue a warning otherwise. + if (converged) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); + } + } + + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ */ \ No newline at end of file diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index 2e9a016f5..fbdd13782 100644 --- a/src/storage/BitVector.h +++ b/src/storage/BitVector.h @@ -149,6 +149,26 @@ public: std::copy(bv.bucketArray, bv.bucketArray + this->bucketCount, this->bucketArray); } + /*! + * Copy Constructor. Performs a deep copy of the bits in the given bit vector that are given by the filter. + * @param bv A reference to the bit vector to be copied. + * @param filter The filter to apply for copying. + */ + BitVector(BitVector const& other, BitVector const& filter) : bitCount(filter.getNumberOfSetBits()), endIterator(*this, bitCount, bitCount, false), truncateMask((1ll << (bitCount & mod64mask)) - 1ll) { + // Determine the number of buckets we need for the given bit count and create bucket array. + bucketCount = bitCount >> 6; + if ((bitCount & mod64mask) != 0) { + ++bucketCount; + } + this->bucketArray = new uint64_t[bucketCount](); + + // Now copy over all bits given by the filter. + uint_fast64_t nextPosition = 0; + for (auto position : filter) { + this->set(nextPosition, other.get(position)); + } + } + /*! * Move constructor. Move constructs the bit vector from the given bit vector. * diff --git a/src/storm.cpp b/src/storm.cpp index 7785122fc..f45f274fc 100644 --- a/src/storm.cpp +++ b/src/storm.cpp @@ -23,8 +23,6 @@ #include "src/storage/SparseMatrix.h" #include "src/models/AtomicPropositionsLabeling.h" #include "src/modelchecker/prctl/EigenDtmcPrctlModelChecker.h" -#include "src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h" -#include "src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h" #include "src/parser/AutoParser.h" #include "src/parser/PrctlParser.h" //#include "src/solver/GraphAnalyzer.h" @@ -144,7 +142,6 @@ void printHeader(const int argc, const char* argv[]) { bool parseOptions(const int argc, const char* argv[]) { storm::settings::Settings* s = nullptr; try { - storm::settings::Settings::registerModule>(); s = storm::settings::newInstance(argc, argv, nullptr); } catch (storm::exceptions::InvalidSettingsException& e) { std::cout << "Could not recover from settings error: " << e.what() << "." << std::endl; @@ -202,7 +199,7 @@ storm::modelchecker::prctl::AbstractModelChecker* createPrctlModelChecke // Create the appropriate model checker. storm::settings::Settings* s = storm::settings::instance(); if (s->getString("matrixlib") == "gmm++") { - return new storm::modelchecker::prctl::GmmxxDtmcPrctlModelChecker(dtmc); + return new storm::modelchecker::prctl::SparseDtmcPrctlModelChecker(dtmc, new storm::solver::AbstractLinearEquationSolver()); } // The control flow should never reach this point, as there is a default setting for matrixlib. @@ -221,9 +218,9 @@ storm::modelchecker::prctl::AbstractModelChecker* createPrctlModelChecke // Create the appropriate model checker. storm::settings::Settings* s = storm::settings::instance(); if (s->getString("matrixlib") == "gmm++") { - return new storm::modelchecker::prctl::GmmxxMdpPrctlModelChecker(mdp); + return new storm::modelchecker::prctl::SparseMdpPrctlModelChecker(mdp, new storm::solver::GmmxxNondeterministicEquationSolver()); } else if (s->getString("matrixlib") == "native") { - return new storm::modelchecker::prctl::SparseMdpPrctlModelChecker(mdp); + return new storm::modelchecker::prctl::SparseMdpPrctlModelChecker(mdp, new storm::solver::AbstractNondeterministicEquationSolver()); } // The control flow should never reach this point, as there is a default setting for matrixlib. diff --git a/src/utility/Settings.cpp b/src/utility/Settings.cpp index d31eea7a2..55e802811 100644 --- a/src/utility/Settings.cpp +++ b/src/utility/Settings.cpp @@ -100,20 +100,20 @@ Settings::Settings(int const argc, char const * const argv[], char const * const } LOG4CPLUS_DEBUG(logger, "Finished loading config."); } - catch (bpo::reading_file e) { + catch (bpo::reading_file const& e) { std::cerr << "Could not read config file " << filename << std::endl; LOG4CPLUS_ERROR(logger, "Could not read config file"); } - catch (bpo::required_option e) { + catch (bpo::required_option const& e) { throw storm::exceptions::InvalidSettingsException() << "Required option missing"; } - catch (bpo::validation_error e) { + catch (bpo::validation_error const& e) { throw storm::exceptions::InvalidSettingsException() << "Validation failed: " << e.what(); } - catch (bpo::invalid_command_line_syntax e) { + catch (bpo::invalid_command_line_syntax const& e) { throw storm::exceptions::InvalidSettingsException() << e.what(); } - catch (bpo::error e) { + catch (bpo::error const& e) { throw storm::exceptions::InvalidSettingsException() << e.what(); } } diff --git a/src/utility/Settings.h b/src/utility/Settings.h index bdf8fef63..b3af50138 100644 --- a/src/utility/Settings.h +++ b/src/utility/Settings.h @@ -125,17 +125,17 @@ namespace settings { * @endcode */ template - static void registerModule() { + static void registerSolver() { // Get trigger values. - std::pair< std::string, std::string > trigger = T::getOptionTrigger(); + std::string const& name = T::getName(); // Build description name. std::stringstream str; - str << "Options for " << T::getModuleName() << " (" << trigger.first << " = " << trigger.second << ")"; + str << "Options for " << name << ":"; std::shared_ptr desc = std::shared_ptr(new bpo::options_description(str.str())); // Put options into description. T::putOptions(desc.get()); // Store module. - Settings::modules[ trigger ] = desc; + Settings::modules[name] = desc; } friend std::ostream& help(std::ostream& os); @@ -185,7 +185,7 @@ namespace settings { /*! * @brief Contains option descriptions for all modules. */ - static std::map< std::pair< std::string, std::string >, std::shared_ptr > modules; + static std::map> modules; /*! * @brief option mapping. diff --git a/src/utility/vector.h b/src/utility/vector.h index 9b973de04..52e16c337 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -100,6 +100,23 @@ void selectVectorValues(std::vector& vector, storm::storage::BitVector const& } } +/*! + * Selects one element out of each row group and writes it to the target vector. + * + * @param vector The target vector to which the values are written. + * @param rowGroupToRowIndexMapping A mapping from row group indices to an offset that specifies which of the values to + * take from the row group. + * @param rowGrouping A vector that specifies the begin and end of each group of elements in the values vector. + * @param values The vector from which to select the values. + */ +template +void selectVectorValues(std::vector& vector, std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGrouping, std::vector const& values) { + uint_fast64_t oldPosition = 0; + for (uint_fast64_t i = 0; i < vector.size(); ++i) { + vector[i] = values[rowGrouping[i] + rowGroupToRowIndexMapping[i]]; + } +} + /*! * Selects values from a vector at the specified positions and writes them into another vector as often as given by * the size of the corresponding group of elements.