From b42aa5f473dd33c9b75b4c30d40eb7e91aba5538 Mon Sep 17 00:00:00 2001 From: TimQu <tim.quatmann@cs.rwth-aachen.de> Date: Thu, 30 Nov 2017 17:30:24 +0100 Subject: [PATCH] initial implementation for quick and sound vi for DTMCs --- .../modules/NativeEquationSolverSettings.cpp | 4 +- .../solver/NativeLinearEquationSolver.cpp | 99 ++++++++++++++++++- src/storm/solver/NativeLinearEquationSolver.h | 1 + src/storm/solver/SolverSelectionOptions.cpp | 2 + src/storm/solver/SolverSelectionOptions.h | 2 +- src/storm/utility/constants.cpp | 6 +- src/storm/utility/constants.h | 2 + .../DtmcPrctlModelCheckerTest.cpp | 19 ++++ .../storm/solver/LinearEquationSolverTest.cpp | 16 +++ 9 files changed, 146 insertions(+), 5 deletions(-) diff --git a/src/storm/settings/modules/NativeEquationSolverSettings.cpp b/src/storm/settings/modules/NativeEquationSolverSettings.cpp index 664818b97..628e0bd75 100644 --- a/src/storm/settings/modules/NativeEquationSolverSettings.cpp +++ b/src/storm/settings/modules/NativeEquationSolverSettings.cpp @@ -25,7 +25,7 @@ namespace storm { const std::string NativeEquationSolverSettings::powerMethodMultiplicationStyleOptionName = "powmult"; NativeEquationSolverSettings::NativeEquationSolverSettings() : ModuleSettings(moduleName) { - std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae", "power", "ratsearch" }; + std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae", "power", "ratsearch", "qpower" }; this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(methods)).setDefaultValueString("jacobi").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); @@ -63,6 +63,8 @@ namespace storm { return storm::solver::NativeLinearEquationSolverMethod::Power; } else if (linearEquationSystemTechniqueAsString == "ratsearch") { return storm::solver::NativeLinearEquationSolverMethod::RationalSearch; + } else if (linearEquationSystemTechniqueAsString == "qpower") { + return storm::solver::NativeLinearEquationSolverMethod::QuickPower; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index c9974a700..d2dc18e68 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -555,6 +555,99 @@ namespace storm { return converged; } + + + template<typename ValueType> + bool NativeLinearEquationSolver<ValueType>::solveEquationsQuickSoundPower(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const { + STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (QuickPower)"); + // Prepare the solution vectors. + if (!this->cachedRowVector) { + this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount(), storm::utility::zero<ValueType>()); + } else { + this->cachedRowVector->assign(getMatrixRowCount(), storm::utility::zero<ValueType>()); + } + if (!this->cachedRowVector2) { + this->cachedRowVector2 = std::make_unique<std::vector<ValueType>>(getMatrixRowCount(), storm::utility::one<ValueType>()); + } else { + this->cachedRowVector2->assign(getMatrixRowCount(), storm::utility::one<ValueType>()); + } + + std::vector<ValueType>* stepBoundedValues = this->cachedRowVector.get(); + std::vector<ValueType>* stepBoundedStayProbabilities = this->cachedRowVector2.get(); + + std::vector<ValueType>* tmp = &x; + + bool converged = false; + bool terminate = false; + uint64_t iterations = 0; + bool doConvergenceCheck = true; + bool useDiffs = this->hasRelevantValues(); + ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision()); + ValueType lowerValueBound, upperValueBound; + bool relative = env.solver().native().getRelativeTerminationCriterion(); + if (!relative) { + precision *= storm::utility::convertNumber<ValueType>(2.0); + } + uint64_t maxIter = env.solver().native().getMaximalNumberOfIterations(); + this->startMeasureProgress(); + auto firstProb1EntryIt = stepBoundedStayProbabilities->begin(); + while (!converged && !terminate && iterations < maxIter) { + this->multiplier.multAdd(*this->A, *stepBoundedValues, &b, *tmp); + std::swap(tmp, stepBoundedValues); + this->multiplier.multAdd(*this->A, *stepBoundedStayProbabilities, nullptr, *tmp); + std::swap(tmp, stepBoundedStayProbabilities); + for (; firstProb1EntryIt != stepBoundedStayProbabilities->end(); ++firstProb1EntryIt) { + static_assert(NumberTraits<ValueType>::IsExact || std::is_same<ValueType, double>::value, "Considered ValueType not handled."); + if (NumberTraits<ValueType>::IsExact) { + if (storm::utility::isOne(*firstProb1EntryIt)) { + break; + } + } else { + if (storm::utility::isAlmostOne(storm::utility::convertNumber<double>(*firstProb1EntryIt))) { + break; + } + } + } + if (firstProb1EntryIt == stepBoundedStayProbabilities->end()) { + auto valIt = stepBoundedValues->begin(); + auto valIte = stepBoundedValues->end(); + auto probIt = stepBoundedStayProbabilities->begin(); + STORM_LOG_ASSERT(!storm::utility::isOne(*probIt), "Did not expect staying-probability 1 at this point."); + lowerValueBound = *valIt / (storm::utility::one<ValueType>() - *probIt); + upperValueBound = lowerValueBound; + ValueType largestStayProb = *probIt; + for (; valIt != valIte; ++valIt, ++probIt) { + ValueType currentBound = *valIt / (storm::utility::one<ValueType>() - *probIt); + lowerValueBound = std::min(lowerValueBound, currentBound); + upperValueBound = std::max(upperValueBound, currentBound); + largestStayProb = std::max(largestStayProb, *probIt); + } + STORM_LOG_ASSERT(!relative, "Relative termination criterion not implemented currently."); + converged = largestStayProb * (upperValueBound - lowerValueBound) < precision; + } + + // Potentially show progress. + this->showProgressIterative(iterations); + + // Set up next iteration. + ++iterations; + + } + + // Finally set up the solution vector + ValueType meanBound = (upperValueBound - lowerValueBound) / storm::utility::convertNumber<ValueType>(2.0); + storm::utility::vector::applyPointwise(*stepBoundedValues, *stepBoundedStayProbabilities, x, [&meanBound] (ValueType const& v, ValueType const& p) { return v + p * meanBound; }); + + if (!this->isCachingEnabled()) { + clearCache(); + } + + this->logIterations(converged, terminate, iterations); + + return converged; + + } + template<typename ValueType> bool NativeLinearEquationSolver<ValueType>::solveEquationsRationalSearch(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const { return solveEquationsRationalSearchHelper<double>(env, x, b); @@ -828,7 +921,7 @@ namespace storm { } else { STORM_LOG_WARN("The selected solution method does not guarantee exact results."); } - } else if (env.solver().isForceSoundness() && method != NativeLinearEquationSolverMethod::Power && method != NativeLinearEquationSolverMethod::RationalSearch) { + } else if (env.solver().isForceSoundness() && method != NativeLinearEquationSolverMethod::Power && method != NativeLinearEquationSolverMethod::RationalSearch && method != NativeLinearEquationSolverMethod::QuickPower) { if (env.solver().native().isMethodSetFromDefault()) { method = NativeLinearEquationSolverMethod::Power; STORM_LOG_INFO("Selecting '" + toString(method) + "' as the solution technique to guarantee sound results. If you want to override this, please explicitly specify a different method."); @@ -857,6 +950,8 @@ namespace storm { } else { return this->solveEquationsPower(env, x, b); } + case NativeLinearEquationSolverMethod::QuickPower: + return this->solveEquationsQuickSoundPower(env, x, b); case NativeLinearEquationSolverMethod::RationalSearch: return this->solveEquationsRationalSearch(env, x, b); } @@ -921,7 +1016,7 @@ namespace storm { template<typename ValueType> LinearEquationSolverProblemFormat NativeLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const { auto method = getMethod(env, storm::NumberTraits<ValueType>::IsExact); - if (method == NativeLinearEquationSolverMethod::Power || method == NativeLinearEquationSolverMethod::RationalSearch) { + if (method == NativeLinearEquationSolverMethod::Power || method == NativeLinearEquationSolverMethod::RationalSearch || method == NativeLinearEquationSolverMethod::QuickPower) { return LinearEquationSolverProblemFormat::FixedPointSystem; } else { return LinearEquationSolverProblemFormat::EquationSystem; diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index 136c8f68f..34ff42bf7 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -71,6 +71,7 @@ namespace storm { virtual bool solveEquationsWalkerChae(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; virtual bool solveEquationsPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; virtual bool solveEquationsSoundPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; + virtual bool solveEquationsQuickSoundPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; virtual bool solveEquationsRationalSearch(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; template<typename RationalType, typename ImpreciseType> diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp index 504f3e2a7..0393e6912 100644 --- a/src/storm/solver/SolverSelectionOptions.cpp +++ b/src/storm/solver/SolverSelectionOptions.cpp @@ -78,6 +78,8 @@ namespace storm { return "Power"; case NativeLinearEquationSolverMethod::RationalSearch: return "RationalSearch"; + case NativeLinearEquationSolverMethod::QuickPower: + return "QuickPower"; } return "invalid"; } diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h index ed4605ecb..e5c7a55eb 100644 --- a/src/storm/solver/SolverSelectionOptions.h +++ b/src/storm/solver/SolverSelectionOptions.h @@ -14,7 +14,7 @@ namespace storm { ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination) ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat) - ExtendEnumsWithSelectionField(NativeLinearEquationSolverMethod, Jacobi, GaussSeidel, SOR, WalkerChae, Power, RationalSearch) + ExtendEnumsWithSelectionField(NativeLinearEquationSolverMethod, Jacobi, GaussSeidel, SOR, WalkerChae, Power, RationalSearch, QuickPower) ExtendEnumsWithSelectionField(GmmxxLinearEquationSolverMethod, Bicgstab, Qmr, Gmres) ExtendEnumsWithSelectionField(GmmxxLinearEquationSolverPreconditioner, Ilu, Diagonal, None) ExtendEnumsWithSelectionField(EigenLinearEquationSolverMethod, SparseLU, Bicgstab, DGmres, Gmres) diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index cd6556fed..2a968422f 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -48,7 +48,11 @@ namespace storm { bool isAlmostZero(double const& a) { return a < 1e-12 && a > -1e-12; } - + + bool isAlmostOne(double const& a) { + return a < (1.0 + 1e-12) && a > (1.0 - 1e-12); + } + template<typename ValueType> bool isConstant(ValueType const&) { return true; diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index 5eb1d1fc5..360128f00 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -44,6 +44,8 @@ namespace storm { bool isZero(ValueType const& a); bool isAlmostZero(double const& a); + + bool isAlmostOne(double const& a); template<typename ValueType> bool isConstant(ValueType const& a); diff --git a/src/test/storm/modelchecker/DtmcPrctlModelCheckerTest.cpp b/src/test/storm/modelchecker/DtmcPrctlModelCheckerTest.cpp index d1dc5a369..a19b04c9c 100644 --- a/src/test/storm/modelchecker/DtmcPrctlModelCheckerTest.cpp +++ b/src/test/storm/modelchecker/DtmcPrctlModelCheckerTest.cpp @@ -222,6 +222,24 @@ namespace { } }; + class SparseNativeQuickSoundPowerEnvironment { + public: + static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models + static const storm::settings::modules::CoreSettings::Engine engine = storm::settings::modules::CoreSettings::Engine::Sparse; + static const bool isExact = false; + typedef double ValueType; + typedef storm::models::sparse::Dtmc<ValueType> ModelType; + static storm::Environment createEnvironment() { + storm::Environment env; + env.solver().setForceSoundness(true); + env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native); + env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::QuickPower); + env.solver().native().setRelativeTerminationCriterion(false); + env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-6)); + return env; + } + }; + class SparseNativeRationalSearchEnvironment { public: static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models @@ -447,6 +465,7 @@ namespace { SparseNativeSorEnvironment, SparseNativePowerEnvironment, SparseNativeSoundPowerEnvironment, + SparseNativeQuickSoundPowerEnvironment, SparseNativeRationalSearchEnvironment, HybridSylvanGmmxxGmresEnvironment, HybridCuddNativeJacobiEnvironment, diff --git a/src/test/storm/solver/LinearEquationSolverTest.cpp b/src/test/storm/solver/LinearEquationSolverTest.cpp index 696654658..b4d7501eb 100644 --- a/src/test/storm/solver/LinearEquationSolverTest.cpp +++ b/src/test/storm/solver/LinearEquationSolverTest.cpp @@ -38,6 +38,21 @@ namespace { } }; + class NativeDoubleQuickSoundPowerEnvironment { + public: + typedef double ValueType; + static const bool isExact = false; + static storm::Environment createEnvironment() { + storm::Environment env; + env.solver().setForceSoundness(true); + env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native); + env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::QuickPower); + env.solver().native().setRelativeTerminationCriterion(false); + env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber, std::string>("1e-6")); + return env; + } + }; + class NativeDoubleJacobiEnvironment { public: typedef double ValueType; @@ -265,6 +280,7 @@ namespace { typedef ::testing::Types< NativeDoublePowerEnvironment, NativeDoubleSoundPowerEnvironment, + NativeDoubleQuickSoundPowerEnvironment, NativeDoubleJacobiEnvironment, NativeDoubleGaussSeidelEnvironment, NativeDoubleSorEnvironment,