diff --git a/CHANGELOG.md b/CHANGELOG.md index 94c0c14b9..4f9972576 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,7 @@ Version 1.0.x - USE_POPCNT removed in favor of FORCE_POPCNT. The popcnt instruction is used if available due to march=native, unless portable is set. Then, using FORCE_POPCNT enables the use of the SSE 4.2 instruction - Parametric model checking is now handled in a separated library/executable called storm-pars -- Support for long-run average rewwards on Markov automata using the value-iteration based approach by Butkova et al. (TACAS 17) +- Support for long-run average rewards on MDPs and Markov automata using a value-iteration based approach - Support for parsing/building models given in the explicit input format of IMCA. ### Version 1.0.1 (2017/4) diff --git a/resources/3rdparty/cudd-3.0.0/cudd/cuddLCache.c b/resources/3rdparty/cudd-3.0.0/cudd/cuddLCache.c index b439187bf..7983e6428 100644 --- a/resources/3rdparty/cudd-3.0.0/cudd/cuddLCache.c +++ b/resources/3rdparty/cudd-3.0.0/cudd/cuddLCache.c @@ -1343,7 +1343,7 @@ cuddHashTableResize( while (item != NULL) { next = item->next; key = item->key; - posn = ddLCHash2(key[0], key[0], shift); + posn = ddLCHash1(key[0], shift); item->next = buckets[posn]; buckets[posn] = item; item = next; diff --git a/resources/3rdparty/cudd-3.0.0/cudd/cuddPriority.c b/resources/3rdparty/cudd-3.0.0/cudd/cuddPriority.c index dfb530f14..78c0633cd 100644 --- a/resources/3rdparty/cudd-3.0.0/cudd/cuddPriority.c +++ b/resources/3rdparty/cudd-3.0.0/cudd/cuddPriority.c @@ -1,3 +1,4 @@ + /** @file @@ -1314,36 +1315,36 @@ Cudd_bddClosestCube( DdNode *g, int *distance) { - DdNode *res, *acube; - CUDD_VALUE_TYPE rdist; + DdNode *res, *acube = NULL; + CUDD_VALUE_TYPE rdist = DD_PLUS_INF_VAL; + CUDD_VALUE_TYPE epsilon = Cudd_ReadEpsilon(dd); - /* Compute the cube and distance as a single ADD. */ do { + /* Compute the cube and distance as a single ADD. */ + Cudd_SetEpsilon(dd,(CUDD_VALUE_TYPE)0.0); dd->reordered = 0; res = cuddBddClosestCube(dd,f,g,CUDD_CONST_INDEX + 1.0); - } while (dd->reordered == 1); - if (res == NULL) { - if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) { - dd->timeoutHandler(dd, dd->tohArg); + Cudd_SetEpsilon(dd,epsilon); + if (dd->reordered == 0) { + if (res == NULL) { + if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) { + dd->timeoutHandler(dd, dd->tohArg); + } + return(NULL); + } + cuddRef(res); + /* Unpack distance and cube. */ + acube = separateCube(dd, res, &rdist); + Cudd_RecursiveDeref(dd, res); } - return(NULL); - } - cuddRef(res); - - /* Unpack distance and cube. */ - do { - dd->reordered = 0; - acube = separateCube(dd, res, &rdist); } while (dd->reordered == 1); if (acube == NULL) { - Cudd_RecursiveDeref(dd, res); if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) { dd->timeoutHandler(dd, dd->tohArg); } - return(NULL); + return(NULL); } cuddRef(acube); - Cudd_RecursiveDeref(dd, res); /* Convert cube from ADD to BDD. */ do { diff --git a/resources/3rdparty/cudd-3.0.0/cudd/cuddUtil.c b/resources/3rdparty/cudd-3.0.0/cudd/cuddUtil.c index 0271fa2bc..45cb503a2 100644 --- a/resources/3rdparty/cudd-3.0.0/cudd/cuddUtil.c +++ b/resources/3rdparty/cudd-3.0.0/cudd/cuddUtil.c @@ -592,10 +592,7 @@ Cudd_SharingSize( // * to deal with functions that depend on more than 1023, but less // * than 2044 variables and don't have too many minterms. // */ -// printf("vars: %i, min_exp: %.400f, %i\n", nvars, DBL_MIN_EXP, DBL_MIN_EXP); -// printf("1: %f, 2: %i \n", (double)(nvars + DBL_MIN_EXP), (nvars + DBL_MIN_EXP)); // max = pow(2.0,(double)(nvars + DBL_MIN_EXP)); -// printf("max: %f, DD_PLUS_INF_VAL: %.50f\n", max, DD_PLUS_INF_VAL); // if (max >= DD_PLUS_INF_VAL) { // return((double)CUDD_OUT_OF_MEM); // } @@ -615,10 +612,8 @@ Cudd_SharingSize( // /* Minterm count is too large to be scaled back. */ // return(DD_PLUS_INF_VAL); // } else { -// printf("in third case: %.400f\n", res); // /* Undo the scaling. */ // res *= pow(2.0,(double)-DBL_MIN_EXP); -// printf("after rescaling: %f\n", res); // return(res); // } // diff --git a/resources/cmake/find_modules/FindZ3.cmake b/resources/cmake/find_modules/FindZ3.cmake index 39bea2558..cf838874d 100644 --- a/resources/cmake/find_modules/FindZ3.cmake +++ b/resources/cmake/find_modules/FindZ3.cmake @@ -10,7 +10,7 @@ # find include dir by searching for a concrete file, which definitely must be in it find_path(Z3_INCLUDE_DIR NAMES z3++.h - PATHS ENV PATH INCLUDE "/usr/local/include/z3/" "${Z3_ROOT}/include" + PATHS ENV PATH INCLUDE "/usr/include/z3" "/usr/local/include/z3/" "${Z3_ROOT}/include" ) # find library diff --git a/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp b/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp index 155b13dec..3d8d18f8b 100644 --- a/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp +++ b/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp @@ -222,14 +222,11 @@ namespace storm { parameterLifter->specifyRegion(region, dirForParameters); // Set up the solver - auto solver = solverFactory->create(parameterLifter->getMatrix()); - if (storm::NumberTraits::IsExact && dynamic_cast*>(solver.get())) { + if (storm::NumberTraits::IsExact && solverFactory->getMinMaxMethod() == storm::solver::MinMaxMethod::ValueIteration) { STORM_LOG_INFO("Parameter Lifting: Setting solution method for exact MinMaxSolver to policy iteration"); - auto* standardSolver = dynamic_cast*>(solver.get()); - auto settings = standardSolver->getSettings(); - settings.setSolutionMethod(storm::solver::StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration); - standardSolver->setSettings(settings); + solverFactory->setMinMaxMethod(storm::solver::MinMaxMethod::PolicyIteration); } + auto solver = solverFactory->create(parameterLifter->getMatrix()); if (lowerResultBound) solver->setLowerBound(lowerResultBound.get()); if (upperResultBound) solver->setUpperBound(upperResultBound.get()); if (!stepBound) solver->setTrackScheduler(true); diff --git a/src/storm/generator/NextStateGenerator.h b/src/storm/generator/NextStateGenerator.h index f894d84a3..25335dbfc 100644 --- a/src/storm/generator/NextStateGenerator.h +++ b/src/storm/generator/NextStateGenerator.h @@ -44,6 +44,8 @@ namespace storm { */ NextStateGenerator(storm::expressions::ExpressionManager const& expressionManager, NextStateGeneratorOptions const& options); + virtual ~NextStateGenerator() = default; + uint64_t getStateSize() const; virtual ModelType getModelType() const = 0; virtual bool isDeterministicModel() const = 0; diff --git a/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp index 1bd826c7f..7ca93bbd4 100644 --- a/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp +++ b/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp @@ -111,7 +111,7 @@ namespace storm { std::unique_ptr SparseMarkovAutomatonCslModelChecker::computeLongRunAverageRewards(storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { STORM_LOG_THROW(checkTask.isOptimizationDirectionSet(), storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); STORM_LOG_THROW(this->getModel().isClosed(), storm::exceptions::InvalidPropertyException, "Unable to compute long run average rewards in non-closed Markov automaton."); - std::vector result = storm::modelchecker::helper::SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(checkTask.getOptimizationDirection(), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getExitRates(), this->getModel().getMarkovianStates(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getRewardModel(""), *minMaxLinearEquationSolverFactory); + std::vector result = storm::modelchecker::helper::SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(checkTask.getOptimizationDirection(), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getExitRates(), this->getModel().getMarkovianStates(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getUniqueRewardModel(), *minMaxLinearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); } diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 722e4c296..1ca0b3728 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -9,6 +9,7 @@ #include "storm/settings/SettingsManager.h" #include "storm/settings/modules/GeneralSettings.h" +#include "storm/settings/modules/MinMaxEquationSolverSettings.h" #include "storm/utility/macros.h" #include "storm/utility/vector.h" @@ -213,7 +214,7 @@ namespace storm { } template - std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { + std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { uint64_t numberOfStates = transitionMatrix.getRowGroupCount(); @@ -233,12 +234,12 @@ namespace storm { storm::utility::vector::setVectorValues(stateRewards, markovianStates & psiStates, storm::utility::one()); storm::models::sparse::StandardRewardModel rewardModel(std::move(stateRewards)); - return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, rewardModel, minMaxLinearEquationSolverFactory, useLpBasedTechnique); + return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, rewardModel, minMaxLinearEquationSolverFactory); } template - std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { + std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { uint64_t numberOfStates = transitionMatrix.getRowGroupCount(); @@ -267,7 +268,7 @@ namespace storm { } // Compute the LRA value for the current MEC. - lraValuesForEndComponents.push_back(computeLraForMaximalEndComponent(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec, minMaxLinearEquationSolverFactory, useLpBasedTechnique)); + lraValuesForEndComponents.push_back(computeLraForMaximalEndComponent(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec, minMaxLinearEquationSolverFactory)); } // For fast transition rewriting, we build some auxiliary data structures. @@ -398,7 +399,7 @@ namespace storm { } template - ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { + ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { // If the mec only consists of a single state, we compute the LRA value directly if (++mec.begin() == mec.end()) { @@ -413,10 +414,14 @@ namespace storm { return result; } - if (useLpBasedTechnique) { + // Solve MEC with the method specified in the settings + storm::solver::LraMethod method = storm::settings::getModule().getLraMethod(); + if (method == storm::solver::LraMethod::LinearProgramming) { return computeLraForMaximalEndComponentLP(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec); - } else { + } else if (method == storm::solver::LraMethod::ValueIteration) { return computeLraForMaximalEndComponentVI(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec, minMaxLinearEquationSolverFactory); + } else { + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } } @@ -560,7 +565,7 @@ namespace storm { // start the iterations - ValueType precision = storm::utility::convertNumber(storm::settings::getModule().getPrecision()) / uniformizationRate; + ValueType precision = storm::utility::convertNumber(storm::settings::getModule().getPrecision()) / uniformizationRate; std::vector v(aMarkovian.getRowCount(), storm::utility::zero()); std::vector w = v; std::vector x(aProbabilistic.getRowGroupCount(), storm::utility::zero()); @@ -610,15 +615,15 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template void SparseMarkovAutomatonCslHelper::computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, double delta, uint64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); @@ -630,15 +635,15 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template void SparseMarkovAutomatonCslHelper::computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, storm::RationalNumber delta, uint64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index f41c8e287..a90c3a0cf 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -29,10 +29,10 @@ namespace storm { template - static std::vector computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); + static std::vector computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template - static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); + static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template static std::vector computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); @@ -61,7 +61,7 @@ namespace storm { * @return The long-run average of being in a goal state for the given MEC. */ template - static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); + static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template static ValueType computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec); template diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index a40d04f08..6b9026ff0 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -43,7 +43,7 @@ namespace storm { template bool SparseMdpPrctlModelChecker::canHandle(CheckTask const& checkTask) const { storm::logic::Formula const& formula = checkTask.getFormula(); - if(formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(false).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true))) { + if(formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true))) { return true; } else { // Check whether we consider a multi-objective formula @@ -162,6 +162,15 @@ namespace storm { return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } + template + std::unique_ptr SparseMdpPrctlModelChecker::computeLongRunAverageRewards(storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { + STORM_LOG_THROW(checkTask.isOptimizationDirectionSet(), storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); + std::vector result = storm::modelchecker::helper::SparseMdpPrctlHelper::computeLongRunAverageRewards(checkTask.getOptimizationDirection(), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getUniqueRewardModel(), *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + } + + + template std::unique_ptr SparseMdpPrctlModelChecker::checkMultiObjectiveFormula(CheckTask const& checkTask) { return multiobjective::performMultiObjectiveModelChecking(this->getModel(), checkTask.getFormula()); diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 91fd30869..43f607031 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -27,6 +27,7 @@ namespace storm { virtual std::unique_ptr computeInstantaneousRewards(storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr computeReachabilityRewards(storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr computeLongRunAverageProbabilities(CheckTask const& checkTask) override; + virtual std::unique_ptr computeLongRunAverageRewards(storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr checkMultiObjectiveFormula(CheckTask const& checkTask) override; private: diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index b10934fdb..c8dab63c4 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -1,4 +1,7 @@ - #include "storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" +#include "storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" + +#include + #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/modelchecker/hints/ExplicitModelCheckerHint.h" @@ -17,6 +20,9 @@ #include "storm/solver/MinMaxLinearEquationSolver.h" #include "storm/solver/LpSolver.h" + +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/MinMaxEquationSolverSettings.h" #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidPropertyException.h" @@ -503,17 +509,31 @@ namespace storm { template std::vector SparseMdpPrctlHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + // If there are no goal states, we avoid the computation and directly return zero. - uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount(); if (psiStates.empty()) { - return std::vector(numberOfStates, storm::utility::zero()); + return std::vector(transitionMatrix.getRowGroupCount(), storm::utility::zero()); } // Likewise, if all bits are set, we can avoid the computation and set. - if ((~psiStates).empty()) { - return std::vector(numberOfStates, storm::utility::one()); + if (psiStates.full()) { + return std::vector(transitionMatrix.getRowGroupCount(), storm::utility::one()); } + // Reduce long run average probabilities to long run average rewards by + // building a reward model assigning one reward to every psi state + std::vector stateRewards(psiStates.size(), storm::utility::zero()); + storm::utility::vector::setVectorValues(stateRewards, psiStates, storm::utility::one()); + storm::models::sparse::StandardRewardModel rewardModel(std::move(stateRewards)); + return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, rewardModel, minMaxLinearEquationSolverFactory); + } + + template + template + std::vector SparseMdpPrctlHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + uint64_t numberOfStates = transitionMatrix.getRowGroupCount(); + // Start by decomposing the MDP into its MECs. storm::storage::MaximalEndComponentDecomposition mecDecomposition(transitionMatrix, backwardTransitions); @@ -529,7 +549,7 @@ namespace storm { for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; - lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(dir, transitionMatrix, psiStates, mec); + lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(dir, transitionMatrix, rewardModel, mec, minMaxLinearEquationSolverFactory); // Gather information for later use. for (auto const& stateChoicesPair : mec) { @@ -555,7 +575,9 @@ namespace storm { // Finally, we are ready to create the SSP matrix and right-hand side of the SSP. std::vector b; - typename storm::storage::SparseMatrixBuilder sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size()); + uint64_t numberOfSspStates = numberOfStatesNotInMecs + mecDecomposition.size(); + + typename storm::storage::SparseMatrixBuilder sspMatrixBuilder(0, numberOfSspStates, 0, false, true, numberOfSspStates); // If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications). uint_fast64_t currentChoice = 0; @@ -596,10 +618,9 @@ namespace storm { boost::container::flat_set const& choicesInMec = stateChoicesPair.second; for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) { - std::vector auxiliaryStateToProbabilityMap(mecDecomposition.size()); - // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state. if (choicesInMec.find(choice) == choicesInMec.end()) { + std::vector auxiliaryStateToProbabilityMap(mecDecomposition.size()); b.push_back(storm::utility::zero()); for (auto element : transitionMatrix.getRow(choice)) { @@ -631,9 +652,9 @@ namespace storm { } // Finalize the matrix and solve the corresponding system of equations. - storm::storage::SparseMatrix sspMatrix = sspMatrixBuilder.build(currentChoice); + storm::storage::SparseMatrix sspMatrix = sspMatrixBuilder.build(currentChoice, numberOfSspStates, numberOfSspStates); - std::vector sspResult(numberOfStatesNotInMecs + mecDecomposition.size()); + std::vector sspResult(numberOfSspStates); std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(std::move(sspMatrix)); solver->solveEquations(dir, sspResult, b); @@ -652,7 +673,134 @@ namespace storm { } template - ValueType SparseMdpPrctlHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { + template + ValueType SparseMdpPrctlHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + // If the mec only consists of a single state, we compute the LRA value directly + if (++mec.begin() == mec.end()) { + uint64_t state = mec.begin()->first; + auto choiceIt = mec.begin()->second.begin(); + ValueType result = rewardModel.getTotalStateActionReward(state, *choiceIt, transitionMatrix); + for (++choiceIt; choiceIt != mec.begin()->second.end(); ++choiceIt) { + if (storm::solver::minimize(dir)) { + result = std::min(result, rewardModel.getTotalStateActionReward(state, *choiceIt, transitionMatrix)); + } else { + result = std::max(result, rewardModel.getTotalStateActionReward(state, *choiceIt, transitionMatrix)); + } + } + return result; + } + + // Solve MEC with the method specified in the settings + storm::solver::LraMethod method = storm::settings::getModule().getLraMethod(); + if (method == storm::solver::LraMethod::LinearProgramming) { + return computeLraForMaximalEndComponentLP(dir, transitionMatrix, rewardModel, mec); + } else if (method == storm::solver::LraMethod::ValueIteration) { + return computeLraForMaximalEndComponentVI(dir, transitionMatrix, rewardModel, mec, minMaxLinearEquationSolverFactory); + } else { + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); + } + } + + template + template + ValueType SparseMdpPrctlHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + // Initialize data about the mec + storm::storage::BitVector mecStates(transitionMatrix.getRowGroupCount(), false); + storm::storage::BitVector mecChoices(transitionMatrix.getRowCount(), false); + for (auto const& stateChoicesPair : mec) { + mecStates.set(stateChoicesPair.first); + for (auto const& choice : stateChoicesPair.second) { + mecChoices.set(choice); + } + } + + boost::container::flat_map toSubModelStateMapping; + uint64_t currState = 0; + toSubModelStateMapping.reserve(mecStates.getNumberOfSetBits()); + for (auto const& mecState : mecStates) { + toSubModelStateMapping.insert(std::pair(mecState, currState)); + ++currState; + } + + // Get a transition matrix that only considers the states and choices within the MEC + storm::storage::SparseMatrixBuilder mecTransitionBuilder(mecChoices.getNumberOfSetBits(), mecStates.getNumberOfSetBits(), 0, true, true, mecStates.getNumberOfSetBits()); + std::vector choiceRewards; + choiceRewards.reserve(mecChoices.getNumberOfSetBits()); + uint64_t currRow = 0; + ValueType selfLoopProb = storm::utility::convertNumber(0.1); // todo try other values + ValueType scalingFactor = storm::utility::one() - selfLoopProb; + for (auto const& mecState : mecStates) { + mecTransitionBuilder.newRowGroup(currRow); + uint64_t groupStart = transitionMatrix.getRowGroupIndices()[mecState]; + uint64_t groupEnd = transitionMatrix.getRowGroupIndices()[mecState + 1]; + for (uint64_t choice = mecChoices.getNextSetIndex(groupStart); choice < groupEnd; choice = mecChoices.getNextSetIndex(choice + 1)) { + bool insertedDiagElement = false; + for (auto const& entry : transitionMatrix.getRow(choice)) { + uint64_t column = toSubModelStateMapping[entry.getColumn()]; + if (!insertedDiagElement && entry.getColumn() > mecState) { + mecTransitionBuilder.addNextValue(currRow, toSubModelStateMapping[mecState], selfLoopProb); + insertedDiagElement = true; + } + if (!insertedDiagElement && entry.getColumn() == mecState) { + mecTransitionBuilder.addNextValue(currRow, column, selfLoopProb + scalingFactor * entry.getValue()); + insertedDiagElement = true; + } else { + mecTransitionBuilder.addNextValue(currRow, column, scalingFactor * entry.getValue()); + } + } + if (!insertedDiagElement) { + mecTransitionBuilder.addNextValue(currRow, toSubModelStateMapping[mecState], selfLoopProb); + } + + // Compute the rewards obtained for this choice + choiceRewards.push_back(scalingFactor * rewardModel.getTotalStateActionReward(mecState, choice, transitionMatrix)); + + ++currRow; + } + } + auto mecTransitions = mecTransitionBuilder.build(); + STORM_LOG_ASSERT(mecTransitions.isProbabilistic(), "The MEC-Matrix is not probabilistic."); + + // start the iterations + ValueType precision = storm::utility::convertNumber(storm::settings::getModule().getPrecision()); + std::vector x(mecTransitions.getRowGroupCount(), storm::utility::zero()); + std::vector xPrime = x; + auto solver = minMaxLinearEquationSolverFactory.create(std::move(mecTransitions)); + solver->setCachingEnabled(true); + ValueType maxDiff, minDiff; + while (true) { + // Compute the obtained rewards for the next step + solver->repeatedMultiply(dir, x, &choiceRewards, 1); + + // update xPrime and check for convergence + // to avoid large (and numerically unstable) x-values, we substract a reference value. + auto xIt = x.begin(); + auto xPrimeIt = xPrime.begin(); + ValueType refVal = *xIt; + maxDiff = *xIt - *xPrimeIt; + minDiff = maxDiff; + *xIt -= refVal; + *xPrimeIt = *xIt; + for (++xIt, ++xPrimeIt; xIt != x.end(); ++xIt, ++xPrimeIt) { + ValueType diff = *xIt - *xPrimeIt; + maxDiff = std::max(maxDiff, diff); + minDiff = std::min(minDiff, diff); + *xIt -= refVal; + *xPrimeIt = *xIt; + } + + if ((maxDiff - minDiff) < precision) { + break; + } + } + return (maxDiff + minDiff) / (storm::utility::convertNumber(2.0) * scalingFactor); + } + + template + template + ValueType SparseMdpPrctlHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec) { std::shared_ptr solver = storm::utility::solver::getLpSolver("LRA for MEC"); solver->setOptimizationDirection(invert(dir)); @@ -672,14 +820,11 @@ namespace storm { // Now, based on the type of the state, create a suitable constraint. for (auto choice : stateChoicesPair.second) { storm::expressions::Expression constraint = -lambda; - ValueType r = 0; for (auto element : transitionMatrix.getRow(choice)) { constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(storm::utility::convertNumber(element.getValue())); - if (psiStates.get(element.getColumn())) { - r += element.getValue(); - } } + typename RewardModelType::ValueType r = rewardModel.getTotalStateActionReward(state, choice, transitionMatrix); constraint = solver->getConstant(storm::utility::convertNumber(r)) + constraint; if (dir == OptimizationDirection::Minimize) { @@ -803,6 +948,10 @@ namespace storm { template std::vector SparseMdpPrctlHelper::computeCumulativeRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint); template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeReachabilityRewards(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint); + template std::vector SparseMdpPrctlHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template double SparseMdpPrctlHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template double SparseMdpPrctlHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template double SparseMdpPrctlHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); #ifdef STORM_HAVE_CARL template class SparseMdpPrctlHelper; @@ -810,6 +959,11 @@ namespace storm { template std::vector SparseMdpPrctlHelper::computeCumulativeRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint); template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeReachabilityRewards(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint); + template std::vector SparseMdpPrctlHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template storm::RationalNumber SparseMdpPrctlHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template storm::RationalNumber SparseMdpPrctlHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template storm::RationalNumber SparseMdpPrctlHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); + #endif } } diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h index 926fa208c..e82c89cc2 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h @@ -65,12 +65,21 @@ namespace storm { static std::vector computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + + template + static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::unique_ptr computeConditionalProbabilities(OptimizationDirection dir, storm::storage::sparse::state_type initialState, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); private: static MDPSparseModelCheckingHelperReturnType computeReachabilityRewardsHelper(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint()); - - static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec); + + template + static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template + static ValueType computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template + static ValueType computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec); }; diff --git a/src/storm/parser/FormulaParserGrammar.cpp b/src/storm/parser/FormulaParserGrammar.cpp index 32148b3c3..33e27e9c8 100644 --- a/src/storm/parser/FormulaParserGrammar.cpp +++ b/src/storm/parser/FormulaParserGrammar.cpp @@ -20,7 +20,7 @@ namespace storm { // Set the identifier mapping to actually generate expressions. expressionParser.setIdentifierMapping(&identifiers_); - longRunAverageRewardFormula = (qi::lit("LRA") | qi::lit("S"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createLongRunAverageRewardFormula, phoenix::ref(*this))]; + longRunAverageRewardFormula = (qi::lit("LRA") | qi::lit("S") | qi::lit("MP"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createLongRunAverageRewardFormula, phoenix::ref(*this))]; longRunAverageRewardFormula.name("long run average reward formula"); instantaneousRewardFormula = (qi::lit("I=") > expressionParser)[qi::_val = phoenix::bind(&FormulaParserGrammar::createInstantaneousRewardFormula, phoenix::ref(*this), qi::_1)]; diff --git a/src/storm/parser/PrismParser.cpp b/src/storm/parser/PrismParser.cpp index 62ef9eebe..0230fbaca 100644 --- a/src/storm/parser/PrismParser.cpp +++ b/src/storm/parser/PrismParser.cpp @@ -443,23 +443,35 @@ namespace storm { } storm::prism::StateReward PrismParser::createStateReward(storm::expressions::Expression statePredicateExpression, storm::expressions::Expression rewardValueExpression) const { - return storm::prism::StateReward(statePredicateExpression, rewardValueExpression, this->getFilename()); + if (this->secondRun) { + return storm::prism::StateReward(statePredicateExpression, rewardValueExpression, this->getFilename()); + } else { + return storm::prism::StateReward(); + } } storm::prism::StateActionReward PrismParser::createStateActionReward(boost::optional const& actionName, storm::expressions::Expression statePredicateExpression, storm::expressions::Expression rewardValueExpression, GlobalProgramInformation& globalProgramInformation) const { - std::string realActionName = actionName ? actionName.get() : ""; - - auto const& nameIndexPair = globalProgramInformation.actionIndices.find(realActionName); - STORM_LOG_THROW(nameIndexPair != globalProgramInformation.actionIndices.end(), storm::exceptions::WrongFormatException, "Transition reward refers to illegal action '" << realActionName << "'."); - return storm::prism::StateActionReward(nameIndexPair->second, realActionName, statePredicateExpression, rewardValueExpression, this->getFilename()); + if (this->secondRun) { + std::string realActionName = actionName ? actionName.get() : ""; + + auto const& nameIndexPair = globalProgramInformation.actionIndices.find(realActionName); + STORM_LOG_THROW(nameIndexPair != globalProgramInformation.actionIndices.end(), storm::exceptions::WrongFormatException, "Transition reward refers to illegal action '" << realActionName << "'."); + return storm::prism::StateActionReward(nameIndexPair->second, realActionName, statePredicateExpression, rewardValueExpression, this->getFilename()); + } else { + return storm::prism::StateActionReward(); + } } storm::prism::TransitionReward PrismParser::createTransitionReward(boost::optional const& actionName, storm::expressions::Expression sourceStatePredicateExpression, storm::expressions::Expression targetStatePredicateExpression, storm::expressions::Expression rewardValueExpression, GlobalProgramInformation& globalProgramInformation) const { - std::string realActionName = actionName ? actionName.get() : ""; - - auto const& nameIndexPair = globalProgramInformation.actionIndices.find(realActionName); - STORM_LOG_THROW(nameIndexPair != globalProgramInformation.actionIndices.end(), storm::exceptions::WrongFormatException, "Transition reward refers to illegal action '" << realActionName << "'."); - return storm::prism::TransitionReward(nameIndexPair->second, realActionName, sourceStatePredicateExpression, targetStatePredicateExpression, rewardValueExpression, this->getFilename()); + if (this->secondRun) { + std::string realActionName = actionName ? actionName.get() : ""; + + auto const& nameIndexPair = globalProgramInformation.actionIndices.find(realActionName); + STORM_LOG_THROW(nameIndexPair != globalProgramInformation.actionIndices.end(), storm::exceptions::WrongFormatException, "Transition reward refers to illegal action '" << realActionName << "'."); + return storm::prism::TransitionReward(nameIndexPair->second, realActionName, sourceStatePredicateExpression, targetStatePredicateExpression, rewardValueExpression, this->getFilename()); + } else { + return storm::prism::TransitionReward(); + } } storm::prism::Assignment PrismParser::createAssignment(std::string const& variableName, storm::expressions::Expression assignedExpression) const { diff --git a/src/storm/settings/modules/CuddSettings.cpp b/src/storm/settings/modules/CuddSettings.cpp index e680276bf..a10ec1bbc 100644 --- a/src/storm/settings/modules/CuddSettings.cpp +++ b/src/storm/settings/modules/CuddSettings.cpp @@ -20,7 +20,7 @@ namespace storm { const std::string CuddSettings::reorderOptionName = "reorder"; CuddSettings::CuddSettings() : ModuleSettings(moduleName) { - this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, true, "Sets the precision used by Cudd.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision up to which to constants are considered to be different.").setDefaultValueDouble(1e-15).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, true, "Sets the precision used by Cudd.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision up to which to constants are considered to be different.").setDefaultValueDouble(1e-15).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0, 1.0)).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, maximalMemoryOptionName, true, "Sets the upper bound of memory available to Cudd in MB.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("value", "The memory available to Cudd (0 means unlimited).").setDefaultValueUnsignedInteger(4096).build()).build()); diff --git a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp index c140e7520..398eaaeac 100644 --- a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp +++ b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp @@ -17,9 +17,10 @@ namespace storm { const std::string MinMaxEquationSolverSettings::maximalIterationsOptionShortName = "i"; const std::string MinMaxEquationSolverSettings::precisionOptionName = "precision"; const std::string MinMaxEquationSolverSettings::absoluteOptionName = "absolute"; + const std::string MinMaxEquationSolverSettings::lraMethodOptionName = "lramethod"; MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) { - std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "acyclic"}; + std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "linear-programming", "lp", "acyclic"}; this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, false, "Sets which min/max linear equation solving technique is preferred.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(minMaxSolvingTechniques)).setDefaultValueString("vi").build()).build()); @@ -28,6 +29,11 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); + + std::vector lraMethods = {"vi", "value-iteration", "linear-programming", "lp"}; + this->addOption(storm::settings::OptionBuilder(moduleName, lraMethodOptionName, false, "Sets which method is preferred for computing long run averages.") + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(lraMethods)).setDefaultValueString("vi").build()).build()); + } storm::solver::MinMaxMethod MinMaxEquationSolverSettings::getMinMaxEquationSolvingMethod() const { @@ -36,6 +42,8 @@ namespace storm { return storm::solver::MinMaxMethod::ValueIteration; } else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") { return storm::solver::MinMaxMethod::PolicyIteration; + } else if (minMaxEquationSolvingTechnique == "linear-programming" || minMaxEquationSolvingTechnique == "lp") { + return storm::solver::MinMaxMethod::LinearProgramming; } else if (minMaxEquationSolvingTechnique == "acyclic") { return storm::solver::MinMaxMethod::Acyclic; } @@ -70,6 +78,17 @@ namespace storm { return this->getOption(absoluteOptionName).getHasOptionBeenSet() ? MinMaxEquationSolverSettings::ConvergenceCriterion::Absolute : MinMaxEquationSolverSettings::ConvergenceCriterion::Relative; } + storm::solver::LraMethod MinMaxEquationSolverSettings::getLraMethod() const { + std::string lraMethodString = this->getOption(lraMethodOptionName).getArgumentByName("name").getValueAsString(); + if (lraMethodString == "value-iteration" || lraMethodString == "vi") { + return storm::solver::LraMethod::ValueIteration; + } else if (lraMethodString == "linear-programming" || lraMethodString == "lp") { + return storm::solver::LraMethod::LinearProgramming; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique '" << lraMethodString << "'."); + } + + } } } diff --git a/src/storm/settings/modules/MinMaxEquationSolverSettings.h b/src/storm/settings/modules/MinMaxEquationSolverSettings.h index 6129b731f..35bee18ac 100644 --- a/src/storm/settings/modules/MinMaxEquationSolverSettings.h +++ b/src/storm/settings/modules/MinMaxEquationSolverSettings.h @@ -75,6 +75,13 @@ namespace storm { */ ConvergenceCriterion getConvergenceCriterion() const; + /*! + * Retrieves the selected long run average method. + * + * @return The selected long run average method. + */ + storm::solver::LraMethod getLraMethod() const; + // The name of the module. static const std::string moduleName; @@ -84,6 +91,7 @@ namespace storm { static const std::string maximalIterationsOptionShortName; static const std::string precisionOptionName; static const std::string absoluteOptionName; + static const std::string lraMethodOptionName; }; } diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp new file mode 100644 index 000000000..d03e356ac --- /dev/null +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -0,0 +1,511 @@ +#include "storm/solver/IterativeMinMaxLinearEquationSolver.h" + +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/MinMaxEquationSolverSettings.h" + +#include "storm/utility/vector.h" +#include "storm/utility/macros.h" +#include "storm/exceptions/InvalidSettingsException.h" +#include "storm/exceptions/InvalidStateException.h" + +namespace storm { + namespace solver { + + template + IterativeMinMaxLinearEquationSolverSettings::IterativeMinMaxLinearEquationSolverSettings() { + // Get the settings object to customize linear solving. + storm::settings::modules::MinMaxEquationSolverSettings const& settings = storm::settings::getModule(); + + maximalNumberOfIterations = settings.getMaximalIterationCount(); + precision = storm::utility::convertNumber(settings.getPrecision()); + relative = settings.getConvergenceCriterion() == storm::settings::modules::MinMaxEquationSolverSettings::ConvergenceCriterion::Relative; + + setSolutionMethod(settings.getMinMaxEquationSolvingMethod()); + + } + + template + void IterativeMinMaxLinearEquationSolverSettings::setSolutionMethod(SolutionMethod const& solutionMethod) { + this->solutionMethod = solutionMethod; + } + + template + void IterativeMinMaxLinearEquationSolverSettings::setSolutionMethod(MinMaxMethod const& solutionMethod) { + switch (solutionMethod) { + case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break; + case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break; + case MinMaxMethod::Acyclic: this->solutionMethod = SolutionMethod::Acyclic; break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique for iterative MinMax linear equation solver."); + } + } + + template + void IterativeMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) { + this->maximalNumberOfIterations = maximalNumberOfIterations; + } + + template + void IterativeMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) { + this->relative = value; + } + + template + void IterativeMinMaxLinearEquationSolverSettings::setPrecision(ValueType precision) { + this->precision = precision; + } + + template + typename IterativeMinMaxLinearEquationSolverSettings::SolutionMethod const& IterativeMinMaxLinearEquationSolverSettings::getSolutionMethod() const { + return solutionMethod; + } + + template + uint64_t IterativeMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const { + return maximalNumberOfIterations; + } + + template + ValueType IterativeMinMaxLinearEquationSolverSettings::getPrecision() const { + return precision; + } + + template + bool IterativeMinMaxLinearEquationSolverSettings::getRelativeTerminationCriterion() const { + return relative; + } + + template + IterativeMinMaxLinearEquationSolver::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, IterativeMinMaxLinearEquationSolverSettings const& settings) : StandardMinMaxLinearEquationSolver(A, std::move(linearEquationSolverFactory)), settings(settings) { + // Intentionally left empty. + } + + template + IterativeMinMaxLinearEquationSolver::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, IterativeMinMaxLinearEquationSolverSettings const& settings) : StandardMinMaxLinearEquationSolver(std::move(A), std::move(linearEquationSolverFactory)), settings(settings) { + // Intentionally left empty. + } + + template + bool IterativeMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + switch (this->getSettings().getSolutionMethod()) { + case IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: + return solveEquationsValueIteration(dir, x, b); + case IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: + return solveEquationsPolicyIteration(dir, x, b); + case IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::Acyclic: + return solveEquationsAcyclic(dir, x, b); + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "This solver does not implement the selected solution method"); + } + return false; + } + + template + bool IterativeMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Create the initial scheduler. + std::vector scheduler = this->hasSchedulerHint() ? this->choicesHint.get() : std::vector(this->A.getRowGroupCount()); + + // Get a vector for storing the right-hand side of the inner equation system. + if(!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A.getRowGroupCount()); + } + std::vector& subB = *auxiliaryRowGroupVector; + + // Resolve the nondeterminism according to the current scheduler. + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); + + // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. + auto solver = this->linearEquationSolverFactory->create(std::move(submatrix)); + if (this->lowerBound) { + solver->setLowerBound(this->lowerBound.get()); + } + if (this->upperBound) { + solver->setUpperBound(this->upperBound.get()); + } + solver->setCachingEnabled(true); + + Status status = Status::InProgress; + uint64_t iterations = 0; + do { + // Solve the equation system for the 'DTMC'. + // FIXME: we need to remove the 0- and 1- states to make the solution unique. + // HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying + // within illegal MECs will never strictly improve the value. Is this true? + solver->solveEquations(x, subB); + + // Go through the multiplication result and see whether we can improve any of the choices. + bool schedulerImproved = false; + for (uint_fast64_t group = 0; group < this->A.getRowGroupCount(); ++group) { + uint_fast64_t currentChoice = scheduler[group]; + for (uint_fast64_t choice = this->A.getRowGroupIndices()[group]; choice < this->A.getRowGroupIndices()[group + 1]; ++choice) { + // If the choice is the currently selected one, we can skip it. + if (choice - this->A.getRowGroupIndices()[group] == currentChoice) { + continue; + } + + // Create the value of the choice. + ValueType choiceValue = storm::utility::zero(); + for (auto const& entry : this->A.getRow(choice)) { + choiceValue += entry.getValue() * x[entry.getColumn()]; + } + choiceValue += b[choice]; + + // If the value is strictly better than the solution of the inner system, we need to improve the scheduler. + // TODO: If the underlying solver is not precise, this might run forever (i.e. when a state has two choices where the (exact) values are equal). + // only changing the scheduler if the values are not equal (modulo precision) would make this unsound. + if (valueImproved(dir, x[group], choiceValue)) { + schedulerImproved = true; + scheduler[group] = choice - this->A.getRowGroupIndices()[group]; + x[group] = std::move(choiceValue); + } + } + } + + // If the scheduler did not improve, we are done. + if (!schedulerImproved) { + status = Status::Converged; + } else { + // Update the scheduler and the solver. + submatrix = this->A.selectRowsFromRowGroups(scheduler, true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); + solver->setMatrix(std::move(submatrix)); + } + + // Update environment variables. + ++iterations; + status = updateStatusIfNotConverged(status, x, iterations); + } while (status == Status::InProgress); + + reportStatus(status, iterations); + + // If requested, we store the scheduler for retrieval. + if (this->isTrackSchedulerSet()) { + this->schedulerChoices = std::move(scheduler); + } + + if(!this->isCachingEnabled()) { + clearCache(); + } + + return status == Status::Converged || status == Status::TerminatedEarly; + } + + template + bool IterativeMinMaxLinearEquationSolver::valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const { + if (dir == OptimizationDirection::Minimize) { + return value2 < value1; + } else { + return value2 > value1; + } + } + + template + ValueType IterativeMinMaxLinearEquationSolver::getPrecision() const { + return this->getSettings().getPrecision(); + } + + template + bool IterativeMinMaxLinearEquationSolver::getRelative() const { + return this->getSettings().getRelativeTerminationCriterion(); + } + + template + bool IterativeMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + if(!this->linEqSolverA) { + this->linEqSolverA = this->linearEquationSolverFactory->create(this->A); + this->linEqSolverA->setCachingEnabled(true); + } + + if (!this->auxiliaryRowVector) { + this->auxiliaryRowVector = std::make_unique>(this->A.getRowCount()); + } + std::vector& multiplyResult = *this->auxiliaryRowVector; + + if (!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A.getRowGroupCount()); + } + + if (this->hasSchedulerHint()) { + // Resolve the nondeterminism according to the scheduler hint + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(this->choicesHint.get(), true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(*auxiliaryRowGroupVector, this->choicesHint.get(), this->A.getRowGroupIndices(), b); + + // Solve the resulting equation system. + // Note that the linEqSolver might consider a slightly different interpretation of "equalModuloPrecision". Hence, we iteratively increase its precision. + auto submatrixSolver = this->linearEquationSolverFactory->create(std::move(submatrix)); + submatrixSolver->setCachingEnabled(true); + if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); } + if (this->upperBound) { submatrixSolver->setUpperBound(this->upperBound.get()); } + submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector); + } + + std::vector* newX = auxiliaryRowGroupVector.get(); + + std::vector* currentX = &x; + + // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. + uint64_t iterations = 0; + + Status status = Status::InProgress; + while (status == Status::InProgress) { + // Compute x' = A*x + b. + this->linEqSolverA->multiply(*currentX, &b, multiplyResult); + + // Reduce the vector x' by applying min/max for all non-deterministic choices. + storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, *newX, this->A.getRowGroupIndices()); + + // Determine whether the method converged. + if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) { + status = Status::Converged; + } + + // Update environment variables. + std::swap(currentX, newX); + ++iterations; + status = updateStatusIfNotConverged(status, *currentX, iterations); + } + + reportStatus(status, 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 (currentX == auxiliaryRowGroupVector.get()) { + std::swap(x, *currentX); + } + + // If requested, we store the scheduler for retrieval. + if (this->isTrackSchedulerSet()) { + // Due to a custom termination condition, it may be the case that no iterations are performed. In this + // case we need to compute x'= A*x+b once. + if (iterations==0) { + this->linEqSolverA->multiply(x, &b, multiplyResult); + } + this->schedulerChoices = std::vector(this->A.getRowGroupCount()); + // Reduce the multiplyResult and keep track of the choices made + storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, x, this->A.getRowGroupIndices(), &this->schedulerChoices.get()); + } + + if (!this->isCachingEnabled()) { + clearCache(); + } + + return status == Status::Converged || status == Status::TerminatedEarly; + } + + template + bool IterativeMinMaxLinearEquationSolver::solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + uint64_t numGroups = this->A.getRowGroupCount(); + + // Allocate memory for the scheduler (if required) + if (this->isTrackSchedulerSet()) { + if (this->schedulerChoices) { + this->schedulerChoices->resize(numGroups); + } else { + this->schedulerChoices = std::vector(numGroups); + } + } + + // We now compute a topological sort of the row groups. + // If caching is enabled, it might be possible to obtain this sort from the cache. + if (this->isCachingEnabled()) { + if (rowGroupOrdering) { + for (auto const& group : *rowGroupOrdering) { + computeOptimalValueForRowGroup(group, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[group] : nullptr); + } + return true; + } else { + rowGroupOrdering = std::make_unique>(); + rowGroupOrdering->reserve(numGroups); + } + } + + auto transposedMatrix = this->A.transpose(true); + + // We store the groups that have already been processed, i.e., the groups for which x[group] was already set to the correct value. + storm::storage::BitVector processedGroups(numGroups, false); + // Furthermore, we keep track of all candidate groups for which we still need to check whether this group can be processed now. + // A group can be processed if all successors have already been processed. + // Notice that the BitVector candidates is considered in a reversed way, i.e., group i is a candidate iff candidates.get(numGroups - i - 1) is true. + // This is due to the observation that groups with higher indices usually need to be processed earlier. + storm::storage::BitVector candidates(numGroups, true); + uint64_t candidate = numGroups - 1; + for (uint64_t numCandidates = candidates.size(); numCandidates > 0; --numCandidates) { + candidates.set(numGroups - candidate - 1, false); + + // Check if the candidate row group has an unprocessed successor + bool hasUnprocessedSuccessor = false; + for (auto const& entry : this->A.getRowGroup(candidate)) { + if (!processedGroups.get(entry.getColumn())) { + hasUnprocessedSuccessor = true; + break; + } + } + + uint64_t nextCandidate = numGroups - candidates.getNextSetIndex(numGroups - candidate - 1 + 1) - 1; + + if (!hasUnprocessedSuccessor) { + // This candidate can be processed. + processedGroups.set(candidate); + computeOptimalValueForRowGroup(candidate, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[candidate] : nullptr); + if (this->isCachingEnabled()) { + rowGroupOrdering->push_back(candidate); + } + + // Add new candidates + for (auto const& predecessorEntry : transposedMatrix.getRow(candidate)) { + uint64_t predecessor = predecessorEntry.getColumn(); + if (!candidates.get(numGroups - predecessor - 1)) { + candidates.set(numGroups - predecessor - 1, true); + nextCandidate = std::max(nextCandidate, predecessor); + ++numCandidates; + } + } + } + candidate = nextCandidate; + } + + assert(candidates.empty()); + STORM_LOG_THROW(processedGroups.full(), storm::exceptions::InvalidOperationException, "The MinMax equation system is not acyclic."); + return true; + } + + template + void IterativeMinMaxLinearEquationSolver::computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice) const { + uint64_t row = this->A.getRowGroupIndices()[group]; + uint64_t groupEnd = this->A.getRowGroupIndices()[group + 1]; + assert(row != groupEnd); + + auto bIt = b.begin() + row; + ValueType& xi = x[group]; + xi = this->A.multiplyRowWithVector(row, x) + *bIt; + uint64_t optimalRow = row; + + for (++row, ++bIt; row < groupEnd; ++row, ++bIt) { + ValueType choiceVal = this->A.multiplyRowWithVector(row, x) + *bIt; + if (minimize(dir)) { + if (choiceVal < xi) { + xi = choiceVal; + optimalRow = row; + } + } else { + if (choiceVal > xi) { + xi = choiceVal; + optimalRow = row; + } + } + } + if (choice != nullptr) { + *choice = optimalRow - this->A.getRowGroupIndices()[group]; + } + } + + template + typename IterativeMinMaxLinearEquationSolver::Status IterativeMinMaxLinearEquationSolver::updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const { + if (status != Status::Converged) { + if (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)) { + status = Status::TerminatedEarly; + } else if (iterations >= this->getSettings().getMaximalNumberOfIterations()) { + status = Status::MaximalIterationsExceeded; + } + } + return status; + } + + template + void IterativeMinMaxLinearEquationSolver::reportStatus(Status status, uint64_t iterations) const { + switch (status) { + case Status::Converged: STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations."); break; + case Status::TerminatedEarly: STORM_LOG_INFO("Iterative solver terminated early after " << iterations << " iterations."); break; + case Status::MaximalIterationsExceeded: STORM_LOG_WARN("Iterative solver did not converge after " << iterations << " iterations."); break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Iterative solver terminated unexpectedly."); + } + } + + template + IterativeMinMaxLinearEquationSolverSettings const& IterativeMinMaxLinearEquationSolver::getSettings() const { + return settings; + } + + template + void IterativeMinMaxLinearEquationSolver::setSettings(IterativeMinMaxLinearEquationSolverSettings const& newSettings) { + settings = newSettings; + } + + template + void IterativeMinMaxLinearEquationSolver::clearCache() const { + auxiliaryRowGroupVector.reset(); + rowGroupOrdering.reset(); + StandardMinMaxLinearEquationSolver::clearCache(); + } + + template + IterativeMinMaxLinearEquationSolverFactory::IterativeMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(method, trackScheduler) { + settings.setSolutionMethod(this->getMinMaxMethod()); + } + + template + IterativeMinMaxLinearEquationSolverFactory::IterativeMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(std::move(linearEquationSolverFactory), method, trackScheduler) { + settings.setSolutionMethod(this->getMinMaxMethod()); + } + + template + IterativeMinMaxLinearEquationSolverFactory::IterativeMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(solverType, method, trackScheduler) { + settings.setSolutionMethod(this->getMinMaxMethod()); + } + + template + std::unique_ptr> IterativeMinMaxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { + STORM_LOG_ASSERT(this->linearEquationSolverFactory, "Linear equation solver factory not initialized."); + + std::unique_ptr> result = std::make_unique>(std::move(matrix), this->linearEquationSolverFactory->clone(), settings); + result->setTrackScheduler(this->isTrackSchedulerSet()); + return result; + } + + template + std::unique_ptr> IterativeMinMaxLinearEquationSolverFactory::create(storm::storage::SparseMatrix&& matrix) const { + STORM_LOG_ASSERT(this->linearEquationSolverFactory, "Linear equation solver factory not initialized."); + + std::unique_ptr> result = std::make_unique>(std::move(matrix), this->linearEquationSolverFactory->clone(), settings); + result->setTrackScheduler(this->isTrackSchedulerSet()); + return result; + } + + template + IterativeMinMaxLinearEquationSolverSettings& IterativeMinMaxLinearEquationSolverFactory::getSettings() { + return settings; + } + + template + IterativeMinMaxLinearEquationSolverSettings const& IterativeMinMaxLinearEquationSolverFactory::getSettings() const { + return settings; + } + + template + void IterativeMinMaxLinearEquationSolverFactory::setMinMaxMethod(MinMaxMethodSelection const& newMethod) { + MinMaxLinearEquationSolverFactory::setMinMaxMethod(newMethod); + settings.setSolutionMethod(this->getMinMaxMethod()); + } + + template + void IterativeMinMaxLinearEquationSolverFactory::setMinMaxMethod(MinMaxMethod const& newMethod) { + MinMaxLinearEquationSolverFactory::setMinMaxMethod(newMethod); + settings.setSolutionMethod(this->getMinMaxMethod()); + } + + template class IterativeMinMaxLinearEquationSolverSettings; + template class IterativeMinMaxLinearEquationSolver; + template class IterativeMinMaxLinearEquationSolverFactory; + +#ifdef STORM_HAVE_CARL + template class IterativeMinMaxLinearEquationSolverSettings; + template class IterativeMinMaxLinearEquationSolver; + template class IterativeMinMaxLinearEquationSolverFactory; +#endif + } +} diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h new file mode 100644 index 000000000..3d4f06156 --- /dev/null +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h @@ -0,0 +1,96 @@ +#pragma once + +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/StandardMinMaxLinearEquationSolver.h" + +namespace storm { + namespace solver { + + template + class IterativeMinMaxLinearEquationSolverSettings { + public: + IterativeMinMaxLinearEquationSolverSettings(); + + enum class SolutionMethod { + ValueIteration, PolicyIteration, Acyclic + }; + + void setSolutionMethod(SolutionMethod const& solutionMethod); + void setSolutionMethod(MinMaxMethod const& solutionMethod); + void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations); + void setRelativeTerminationCriterion(bool value); + void setPrecision(ValueType precision); + + SolutionMethod const& getSolutionMethod() const; + uint64_t getMaximalNumberOfIterations() const; + ValueType getPrecision() const; + bool getRelativeTerminationCriterion() const; + + private: + SolutionMethod solutionMethod; + uint64_t maximalNumberOfIterations; + ValueType precision; + bool relative; + }; + + template + class IterativeMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver { + public: + IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, IterativeMinMaxLinearEquationSolverSettings const& settings = IterativeMinMaxLinearEquationSolverSettings()); + IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, IterativeMinMaxLinearEquationSolverSettings const& settings = IterativeMinMaxLinearEquationSolverSettings()); + + virtual bool solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; + + IterativeMinMaxLinearEquationSolverSettings const& getSettings() const; + void setSettings(IterativeMinMaxLinearEquationSolverSettings const& newSettings); + + virtual void clearCache() const override; + + virtual ValueType getPrecision() const override; + virtual bool getRelative() const override; + private: + bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + bool solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + bool solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + + bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; + + void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; + + enum class Status { + Converged, TerminatedEarly, MaximalIterationsExceeded, InProgress + }; + + // possibly cached data + mutable std::unique_ptr> auxiliaryRowGroupVector; // A.rowGroupCount() entries + mutable std::unique_ptr> rowGroupOrdering; // A.rowGroupCount() entries + + Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const; + void reportStatus(Status status, uint64_t iterations) const; + + /// The settings of this solver. + IterativeMinMaxLinearEquationSolverSettings settings; + }; + + template + class IterativeMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory { + public: + IterativeMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); + IterativeMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); + IterativeMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); + + virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; + virtual std::unique_ptr> create(storm::storage::SparseMatrix&& matrix) const override; + + IterativeMinMaxLinearEquationSolverSettings& getSettings(); + IterativeMinMaxLinearEquationSolverSettings const& getSettings() const; + + virtual void setMinMaxMethod(MinMaxMethodSelection const& newMethod) override; + virtual void setMinMaxMethod(MinMaxMethod const& newMethod) override; + + private: + IterativeMinMaxLinearEquationSolverSettings settings; + + }; + } +} diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index 4bdb08f3e..1c6f8dc65 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -3,7 +3,7 @@ #include #include "storm/solver/LinearEquationSolver.h" -#include "storm/solver/StandardMinMaxLinearEquationSolver.h" +#include "storm/solver/IterativeMinMaxLinearEquationSolver.h" #include "storm/solver/TopologicalMinMaxLinearEquationSolver.h" #include "storm/settings/SettingsManager.h" @@ -131,8 +131,8 @@ namespace storm { } template - MinMaxLinearEquationSolverFactory::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) { - // Intentionally left empty. + MinMaxLinearEquationSolverFactory::MinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : trackScheduler(trackScheduler) { + setMinMaxMethod(method); } template @@ -151,7 +151,26 @@ namespace storm { } template - GeneralMinMaxLinearEquationSolverFactory::GeneralMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { + void MinMaxLinearEquationSolverFactory::setMinMaxMethod(MinMaxMethodSelection const& newMethod) { + if (newMethod == MinMaxMethodSelection::FROMSETTINGS) { + setMinMaxMethod(storm::settings::getModule().getMinMaxEquationSolvingMethod()); + } else { + setMinMaxMethod(convert(newMethod)); + } + } + + template + void MinMaxLinearEquationSolverFactory::setMinMaxMethod(MinMaxMethod const& newMethod) { + method = newMethod; + } + + template + MinMaxMethod const& MinMaxLinearEquationSolverFactory::getMinMaxMethod() const { + return method; + } + + template + GeneralMinMaxLinearEquationSolverFactory::GeneralMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : MinMaxLinearEquationSolverFactory(method, trackScheduler) { // Intentionally left empty. } @@ -169,9 +188,11 @@ namespace storm { template std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { std::unique_ptr> result; - auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); + auto method = this->getMinMaxMethod(); if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) { - result = std::make_unique>(std::forward(matrix), std::make_unique>()); + IterativeMinMaxLinearEquationSolverSettings iterativeSolverSettings; + iterativeSolverSettings.setSolutionMethod(method); + result = std::make_unique>(std::forward(matrix), std::make_unique>(), iterativeSolverSettings); } else if (method == MinMaxMethod::Topological) { result = std::make_unique>(std::forward(matrix)); } else { @@ -186,9 +207,11 @@ namespace storm { template std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { std::unique_ptr> result; - auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); + auto method = this->getMinMaxMethod(); STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic, storm::exceptions::InvalidSettingsException, "For this data type only value iteration, policy iteration, and acyclic value iteration are available."); - return std::make_unique>(std::forward(matrix), std::make_unique>()); + IterativeMinMaxLinearEquationSolverSettings iterativeSolverSettings; + iterativeSolverSettings.setSolutionMethod(method); + return std::make_unique>(std::forward(matrix), std::make_unique>(), iterativeSolverSettings); } #endif template class MinMaxLinearEquationSolver; diff --git a/src/storm/solver/MinMaxLinearEquationSolver.h b/src/storm/solver/MinMaxLinearEquationSolver.h index 1b869d165..7b4aa18a1 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.h +++ b/src/storm/solver/MinMaxLinearEquationSolver.h @@ -196,22 +196,28 @@ namespace storm { template class MinMaxLinearEquationSolverFactory { public: - MinMaxLinearEquationSolverFactory(bool trackScheduler = false); + MinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const = 0; virtual std::unique_ptr> create(storm::storage::SparseMatrix&& matrix) const; void setTrackScheduler(bool value); bool isTrackSchedulerSet() const; + + virtual void setMinMaxMethod(MinMaxMethodSelection const& newMethod); + virtual void setMinMaxMethod(MinMaxMethod const& newMethod); + + MinMaxMethod const& getMinMaxMethod() const; private: bool trackScheduler; + MinMaxMethod method; }; template class GeneralMinMaxLinearEquationSolverFactory : public MinMaxLinearEquationSolverFactory { public: - GeneralMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + GeneralMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; virtual std::unique_ptr> create(storm::storage::SparseMatrix&& matrix) const override; diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp index 2a39108ea..5a516c8dd 100644 --- a/src/storm/solver/SolverSelectionOptions.cpp +++ b/src/storm/solver/SolverSelectionOptions.cpp @@ -8,6 +8,8 @@ namespace storm { return "policy"; case MinMaxMethod::ValueIteration: return "value"; + case MinMaxMethod::LinearProgramming: + return "linearprogramming"; case MinMaxMethod::Topological: return "topological"; case MinMaxMethod::Acyclic: @@ -16,6 +18,16 @@ namespace storm { return "invalid"; } + std::string toString(LraMethod m) { + switch(m) { + case LraMethod::LinearProgramming: + return "linearprogramming"; + case LraMethod::ValueIteration: + return "valueiteration"; + } + return "invalid"; + } + std::string toString(LpSolverType t) { switch(t) { case LpSolverType::Gurobi: diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h index 94584c7c3..141bc8a0e 100644 --- a/src/storm/solver/SolverSelectionOptions.h +++ b/src/storm/solver/SolverSelectionOptions.h @@ -6,8 +6,9 @@ namespace storm { namespace solver { - ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, Topological, Acyclic) + ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, Acyclic) ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration) + ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk, Z3) ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination) diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp index 627dcd131..6b0f140ce 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp @@ -1,8 +1,6 @@ #include "storm/solver/StandardMinMaxLinearEquationSolver.h" -#include "storm/settings/SettingsManager.h" -#include "storm/settings/modules/MinMaxEquationSolverSettings.h" - +#include "storm/solver/IterativeMinMaxLinearEquationSolver.h" #include "storm/solver/GmmxxLinearEquationSolver.h" #include "storm/solver/EigenLinearEquationSolver.h" #include "storm/solver/NativeLinearEquationSolver.h" @@ -17,395 +15,18 @@ namespace storm { namespace solver { template - StandardMinMaxLinearEquationSolverSettings::StandardMinMaxLinearEquationSolverSettings() { - // Get the settings object to customize linear solving. - storm::settings::modules::MinMaxEquationSolverSettings const& settings = storm::settings::getModule(); - - maximalNumberOfIterations = settings.getMaximalIterationCount(); - precision = storm::utility::convertNumber(settings.getPrecision()); - relative = settings.getConvergenceCriterion() == storm::settings::modules::MinMaxEquationSolverSettings::ConvergenceCriterion::Relative; - - auto method = settings.getMinMaxEquationSolvingMethod(); - switch (method) { - case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break; - case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break; - case MinMaxMethod::Acyclic: this->solutionMethod = SolutionMethod::Acyclic; break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); - } - } - - template - void StandardMinMaxLinearEquationSolverSettings::setSolutionMethod(SolutionMethod const& solutionMethod) { - this->solutionMethod = solutionMethod; - } - - template - void StandardMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) { - this->maximalNumberOfIterations = maximalNumberOfIterations; - } - - template - void StandardMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) { - this->relative = value; - } - - template - void StandardMinMaxLinearEquationSolverSettings::setPrecision(ValueType precision) { - this->precision = precision; - } - - template - typename StandardMinMaxLinearEquationSolverSettings::SolutionMethod const& StandardMinMaxLinearEquationSolverSettings::getSolutionMethod() const { - return solutionMethod; - } - - template - uint64_t StandardMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const { - return maximalNumberOfIterations; - } - - template - ValueType StandardMinMaxLinearEquationSolverSettings::getPrecision() const { - return precision; - } - - template - bool StandardMinMaxLinearEquationSolverSettings::getRelativeTerminationCriterion() const { - return relative; - } - - template - StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(nullptr), A(A) { + StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory) : linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(nullptr), A(A) { // Intentionally left empty. } template - StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(std::make_unique>(std::move(A))), A(*localA) { + StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory) : linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(std::make_unique>(std::move(A))), A(*localA) { // Intentionally left empty. } - template - bool StandardMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - switch (this->getSettings().getSolutionMethod()) { - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: - return solveEquationsValueIteration(dir, x, b); - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: - return solveEquationsPolicyIteration(dir, x, b); - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::Acyclic: - return solveEquationsAcyclic(dir, x, b); - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "This solver does not implement the selected solution method"); - } - return false; - } - - template - bool StandardMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - // Create the initial scheduler. - std::vector scheduler = this->hasSchedulerHint() ? this->choicesHint.get() : std::vector(this->A.getRowGroupCount()); - - // Get a vector for storing the right-hand side of the inner equation system. - if(!auxiliaryRowGroupVector) { - auxiliaryRowGroupVector = std::make_unique>(this->A.getRowGroupCount()); - } - std::vector& subB = *auxiliaryRowGroupVector; - - // Resolve the nondeterminism according to the current scheduler. - storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); - submatrix.convertToEquationSystem(); - storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); - - // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. - auto solver = linearEquationSolverFactory->create(std::move(submatrix)); - if (this->lowerBound) { - solver->setLowerBound(this->lowerBound.get()); - } - if (this->upperBound) { - solver->setUpperBound(this->upperBound.get()); - } - solver->setCachingEnabled(true); - - Status status = Status::InProgress; - uint64_t iterations = 0; - do { - // Solve the equation system for the 'DTMC'. - // FIXME: we need to remove the 0- and 1- states to make the solution unique. - // HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying - // within illegal MECs will never strictly improve the value. Is this true? - solver->solveEquations(x, subB); - - // Go through the multiplication result and see whether we can improve any of the choices. - bool schedulerImproved = false; - for (uint_fast64_t group = 0; group < this->A.getRowGroupCount(); ++group) { - uint_fast64_t currentChoice = scheduler[group]; - for (uint_fast64_t choice = this->A.getRowGroupIndices()[group]; choice < this->A.getRowGroupIndices()[group + 1]; ++choice) { - // If the choice is the currently selected one, we can skip it. - if (choice - this->A.getRowGroupIndices()[group] == currentChoice) { - continue; - } - - // Create the value of the choice. - ValueType choiceValue = storm::utility::zero(); - for (auto const& entry : this->A.getRow(choice)) { - choiceValue += entry.getValue() * x[entry.getColumn()]; - } - choiceValue += b[choice]; - - // If the value is strictly better than the solution of the inner system, we need to improve the scheduler. - // TODO: If the underlying solver is not precise, this might run forever (i.e. when a state has two choices where the (exact) values are equal). - // only changing the scheduler if the values are not equal (modulo precision) would make this unsound. - if (valueImproved(dir, x[group], choiceValue)) { - schedulerImproved = true; - scheduler[group] = choice - this->A.getRowGroupIndices()[group]; - x[group] = std::move(choiceValue); - } - } - } - - // If the scheduler did not improve, we are done. - if (!schedulerImproved) { - status = Status::Converged; - } else { - // Update the scheduler and the solver. - submatrix = this->A.selectRowsFromRowGroups(scheduler, true); - submatrix.convertToEquationSystem(); - storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); - solver->setMatrix(std::move(submatrix)); - } - - // Update environment variables. - ++iterations; - status = updateStatusIfNotConverged(status, x, iterations); - } while (status == Status::InProgress); - - reportStatus(status, iterations); - - // If requested, we store the scheduler for retrieval. - if (this->isTrackSchedulerSet()) { - this->schedulerChoices = std::move(scheduler); - } - - if(!this->isCachingEnabled()) { - clearCache(); - } - - return status == Status::Converged || status == Status::TerminatedEarly; - } - - template - bool StandardMinMaxLinearEquationSolver::valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const { - if (dir == OptimizationDirection::Minimize) { - return value2 < value1; - } else { - return value2 > value1; - } - } - - template - ValueType StandardMinMaxLinearEquationSolver::getPrecision() const { - return this->getSettings().getPrecision(); - } - - template - bool StandardMinMaxLinearEquationSolver::getRelative() const { - return this->getSettings().getRelativeTerminationCriterion(); - } - - template - bool StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - if(!linEqSolverA) { - linEqSolverA = linearEquationSolverFactory->create(A); - linEqSolverA->setCachingEnabled(true); - } - - if (!auxiliaryRowVector) { - auxiliaryRowVector = std::make_unique>(A.getRowCount()); - } - std::vector& multiplyResult = *auxiliaryRowVector; - - if (!auxiliaryRowGroupVector) { - auxiliaryRowGroupVector = std::make_unique>(A.getRowGroupCount()); - } - - if (this->hasSchedulerHint()) { - // Resolve the nondeterminism according to the scheduler hint - storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(this->choicesHint.get(), true); - submatrix.convertToEquationSystem(); - storm::utility::vector::selectVectorValues(*auxiliaryRowGroupVector, this->choicesHint.get(), this->A.getRowGroupIndices(), b); - - // Solve the resulting equation system. - // Note that the linEqSolver might consider a slightly different interpretation of "equalModuloPrecision". Hence, we iteratively increase its precision. - auto submatrixSolver = linearEquationSolverFactory->create(std::move(submatrix)); - submatrixSolver->setCachingEnabled(true); - if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); } - if (this->upperBound) { submatrixSolver->setUpperBound(this->upperBound.get()); } - submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector); - } - - std::vector* newX = auxiliaryRowGroupVector.get(); - - std::vector* currentX = &x; - - // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. - uint64_t iterations = 0; - - Status status = Status::InProgress; - while (status == Status::InProgress) { - // Compute x' = A*x + b. - linEqSolverA->multiply(*currentX, &b, multiplyResult); - - // Reduce the vector x' by applying min/max for all non-deterministic choices. - storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, *newX, this->A.getRowGroupIndices()); - - // Determine whether the method converged. - if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) { - status = Status::Converged; - } - - // Update environment variables. - std::swap(currentX, newX); - ++iterations; - status = updateStatusIfNotConverged(status, *currentX, iterations); - } - - reportStatus(status, 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 (currentX == auxiliaryRowGroupVector.get()) { - std::swap(x, *currentX); - } - - // If requested, we store the scheduler for retrieval. - if (this->isTrackSchedulerSet()) { - // Due to a custom termination condition, it may be the case that no iterations are performed. In this - // case we need to compute x'= A*x+b once. - if (iterations==0) { - linEqSolverA->multiply(x, &b, multiplyResult); - } - this->schedulerChoices = std::vector(this->A.getRowGroupCount()); - // Reduce the multiplyResult and keep track of the choices made - storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, x, this->A.getRowGroupIndices(), &this->schedulerChoices.get()); - } - - if (!this->isCachingEnabled()) { - clearCache(); - } - - return status == Status::Converged || status == Status::TerminatedEarly; - } - - template - bool StandardMinMaxLinearEquationSolver::solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - uint64_t numGroups = this->A.getRowGroupCount(); - - // Allocate memory for the scheduler (if required) - if (this->isTrackSchedulerSet()) { - if (this->schedulerChoices) { - this->schedulerChoices->resize(numGroups); - } else { - this->schedulerChoices = std::vector(numGroups); - } - } - - // We now compute a topological sort of the row groups. - // If caching is enabled, it might be possible to obtain this sort from the cache. - if (this->isCachingEnabled()) { - if (rowGroupOrdering) { - for (auto const& group : *rowGroupOrdering) { - computeOptimalValueForRowGroup(group, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[group] : nullptr); - } - return true; - } else { - rowGroupOrdering = std::make_unique>(); - rowGroupOrdering->reserve(numGroups); - } - } - - auto transposedMatrix = this->A.transpose(true); - - // We store the groups that have already been processed, i.e., the groups for which x[group] was already set to the correct value. - storm::storage::BitVector processedGroups(numGroups, false); - // Furthermore, we keep track of all candidate groups for which we still need to check whether this group can be processed now. - // A group can be processed if all successors have already been processed. - // Notice that the BitVector candidates is considered in a reversed way, i.e., group i is a candidate iff candidates.get(numGroups - i - 1) is true. - // This is due to the observation that groups with higher indices usually need to be processed earlier. - storm::storage::BitVector candidates(numGroups, true); - uint64_t candidate = numGroups - 1; - for (uint64_t numCandidates = candidates.size(); numCandidates > 0; --numCandidates) { - candidates.set(numGroups - candidate - 1, false); - - // Check if the candidate row group has an unprocessed successor - bool hasUnprocessedSuccessor = false; - for (auto const& entry : this->A.getRowGroup(candidate)) { - if (!processedGroups.get(entry.getColumn())) { - hasUnprocessedSuccessor = true; - break; - } - } - - uint64_t nextCandidate = numGroups - candidates.getNextSetIndex(numGroups - candidate - 1 + 1) - 1; - - if (!hasUnprocessedSuccessor) { - // This candidate can be processed. - processedGroups.set(candidate); - computeOptimalValueForRowGroup(candidate, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[candidate] : nullptr); - if (this->isCachingEnabled()) { - rowGroupOrdering->push_back(candidate); - } - - // Add new candidates - for (auto const& predecessorEntry : transposedMatrix.getRow(candidate)) { - uint64_t predecessor = predecessorEntry.getColumn(); - if (!candidates.get(numGroups - predecessor - 1)) { - candidates.set(numGroups - predecessor - 1, true); - nextCandidate = std::max(nextCandidate, predecessor); - ++numCandidates; - } - } - } - candidate = nextCandidate; - } - - assert(candidates.empty()); - STORM_LOG_THROW(processedGroups.full(), storm::exceptions::InvalidOperationException, "The MinMax equation system is not acyclic."); - return true; - } - - template - void StandardMinMaxLinearEquationSolver::computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice) const { - uint64_t row = this->A.getRowGroupIndices()[group]; - uint64_t groupEnd = this->A.getRowGroupIndices()[group + 1]; - assert(row != groupEnd); - - auto bIt = b.begin() + row; - ValueType& xi = x[group]; - xi = this->A.multiplyRowWithVector(row, x) + *bIt; - uint64_t optimalRow = row; - - for (++row, ++bIt; row < groupEnd; ++row, ++bIt) { - ValueType choiceVal = this->A.multiplyRowWithVector(row, x) + *bIt; - if (minimize(dir)) { - if (choiceVal < xi) { - xi = choiceVal; - optimalRow = row; - } - } else { - if (choiceVal > xi) { - xi = choiceVal; - optimalRow = row; - } - } - } - if (choice != nullptr) { - *choice = optimalRow - this->A.getRowGroupIndices()[group]; - } - } - template void StandardMinMaxLinearEquationSolver::repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector const* b, uint_fast64_t n) const { - if(!linEqSolverA) { + if (!linEqSolverA) { linEqSolverA = linearEquationSolverFactory->create(A); linEqSolverA->setCachingEnabled(true); } @@ -423,65 +44,30 @@ namespace storm { storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, x, this->A.getRowGroupIndices()); } - if(!this->isCachingEnabled()) { + if (!this->isCachingEnabled()) { clearCache(); } } - template - typename StandardMinMaxLinearEquationSolver::Status StandardMinMaxLinearEquationSolver::updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const { - if (status != Status::Converged) { - if (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)) { - status = Status::TerminatedEarly; - } else if (iterations >= this->getSettings().getMaximalNumberOfIterations()) { - status = Status::MaximalIterationsExceeded; - } - } - return status; - } - - template - void StandardMinMaxLinearEquationSolver::reportStatus(Status status, uint64_t iterations) const { - switch (status) { - case Status::Converged: STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations."); break; - case Status::TerminatedEarly: STORM_LOG_INFO("Iterative solver terminated early after " << iterations << " iterations."); break; - case Status::MaximalIterationsExceeded: STORM_LOG_WARN("Iterative solver did not converge after " << iterations << " iterations."); break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Iterative solver terminated unexpectedly."); - } - } - - template - StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolver::getSettings() const { - return settings; - } - - template - void StandardMinMaxLinearEquationSolver::setSettings(StandardMinMaxLinearEquationSolverSettings const& newSettings) { - settings = newSettings; - } - template void StandardMinMaxLinearEquationSolver::clearCache() const { linEqSolverA.reset(); auxiliaryRowVector.reset(); - auxiliaryRowGroupVector.reset(); - rowGroupOrdering.reset(); MinMaxLinearEquationSolver::clearCache(); } template - StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler), linearEquationSolverFactory(nullptr) { + StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : MinMaxLinearEquationSolverFactory(method, trackScheduler), linearEquationSolverFactory(std::make_unique>()) { // Intentionally left empty. } template - StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { + StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, MinMaxMethodSelection const& method, bool trackScheduler) : MinMaxLinearEquationSolverFactory(method, trackScheduler), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { // Intentionally left empty. } template - StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { + StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, MinMaxMethodSelection const& method, bool trackScheduler) : MinMaxLinearEquationSolverFactory(method, trackScheduler) { switch (solverType) { case EquationSolverType::Gmmxx: linearEquationSolverFactory = std::make_unique>(); break; case EquationSolverType::Eigen: linearEquationSolverFactory = std::make_unique>(); break; @@ -492,7 +78,7 @@ namespace storm { #ifdef STORM_HAVE_CARL template<> - StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { + StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, MinMaxMethodSelection const& method, bool trackScheduler) : MinMaxLinearEquationSolverFactory(method, trackScheduler) { switch (solverType) { case EquationSolverType::Eigen: linearEquationSolverFactory = std::make_unique>(); break; case EquationSolverType::Elimination: linearEquationSolverFactory = std::make_unique>(); break; @@ -504,63 +90,58 @@ namespace storm { template std::unique_ptr> StandardMinMaxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { + STORM_LOG_ASSERT(linearEquationSolverFactory, "Linear equation solver factory not initialized."); + std::unique_ptr> result; - if (linearEquationSolverFactory) { - result = std::make_unique>(matrix, linearEquationSolverFactory->clone(), settings); + auto method = this->getMinMaxMethod(); + if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) { + IterativeMinMaxLinearEquationSolverSettings iterativeSolverSettings; + iterativeSolverSettings.setSolutionMethod(method); + result = std::make_unique>(matrix, linearEquationSolverFactory->clone(), iterativeSolverSettings); } else { - result = std::make_unique>(matrix, std::make_unique>(), settings); - } - if (this->isTrackSchedulerSet()) { - result->setTrackScheduler(true); + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } + result->setTrackScheduler(this->isTrackSchedulerSet()); return result; } template std::unique_ptr> StandardMinMaxLinearEquationSolverFactory::create(storm::storage::SparseMatrix&& matrix) const { + STORM_LOG_ASSERT(linearEquationSolverFactory, "Linear equation solver factory not initialized."); + std::unique_ptr> result; - if (linearEquationSolverFactory) { - result = std::make_unique>(std::move(matrix), linearEquationSolverFactory->clone(), settings); + auto method = this->getMinMaxMethod(); + if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) { + IterativeMinMaxLinearEquationSolverSettings iterativeSolverSettings; + iterativeSolverSettings.setSolutionMethod(method); + result = std::make_unique>(std::move(matrix), linearEquationSolverFactory->clone(), iterativeSolverSettings); } else { - result = std::make_unique>(std::move(matrix), std::make_unique>(), settings); - } - if (this->isTrackSchedulerSet()) { - result->setTrackScheduler(true); + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } + result->setTrackScheduler(this->isTrackSchedulerSet()); return result; } template - StandardMinMaxLinearEquationSolverSettings& StandardMinMaxLinearEquationSolverFactory::getSettings() { - return settings; - } - - template - StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolverFactory::getSettings() const { - return settings; - } - - template - GmmxxMinMaxLinearEquationSolverFactory::GmmxxMinMaxLinearEquationSolverFactory(bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Gmmxx, trackScheduler) { + GmmxxMinMaxLinearEquationSolverFactory::GmmxxMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Gmmxx, method, trackScheduler) { // Intentionally left empty. } template - EigenMinMaxLinearEquationSolverFactory::EigenMinMaxLinearEquationSolverFactory(bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Eigen, trackScheduler) { + EigenMinMaxLinearEquationSolverFactory::EigenMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Eigen, method, trackScheduler) { // Intentionally left empty. } template - NativeMinMaxLinearEquationSolverFactory::NativeMinMaxLinearEquationSolverFactory(bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Native, trackScheduler) { + NativeMinMaxLinearEquationSolverFactory::NativeMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Native, method, trackScheduler) { // Intentionally left empty. } template - EliminationMinMaxLinearEquationSolverFactory::EliminationMinMaxLinearEquationSolverFactory(bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Elimination, trackScheduler) { + EliminationMinMaxLinearEquationSolverFactory::EliminationMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory(EquationSolverType::Elimination, method, trackScheduler) { // Intentionally left empty. } - template class StandardMinMaxLinearEquationSolverSettings; template class StandardMinMaxLinearEquationSolver; template class StandardMinMaxLinearEquationSolverFactory; template class GmmxxMinMaxLinearEquationSolverFactory; @@ -569,7 +150,6 @@ namespace storm { template class EliminationMinMaxLinearEquationSolverFactory; #ifdef STORM_HAVE_CARL - template class StandardMinMaxLinearEquationSolverSettings; template class StandardMinMaxLinearEquationSolver; template class StandardMinMaxLinearEquationSolverFactory; template class EigenMinMaxLinearEquationSolverFactory; diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.h b/src/storm/solver/StandardMinMaxLinearEquationSolver.h index 050ce0bde..1305f3e26 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.h @@ -6,72 +6,23 @@ namespace storm { namespace solver { - template - class StandardMinMaxLinearEquationSolverSettings { - public: - StandardMinMaxLinearEquationSolverSettings(); - - enum class SolutionMethod { - ValueIteration, PolicyIteration, Acyclic - }; - - void setSolutionMethod(SolutionMethod const& solutionMethod); - void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations); - void setRelativeTerminationCriterion(bool value); - void setPrecision(ValueType precision); - - SolutionMethod const& getSolutionMethod() const; - uint64_t getMaximalNumberOfIterations() const; - ValueType getPrecision() const; - bool getRelativeTerminationCriterion() const; - - private: - SolutionMethod solutionMethod; - uint64_t maximalNumberOfIterations; - ValueType precision; - bool relative; - }; - template class StandardMinMaxLinearEquationSolver : public MinMaxLinearEquationSolver { public: - StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); - StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); + StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory); + StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory); + + virtual ~StandardMinMaxLinearEquationSolver() = default; - virtual bool solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; virtual void repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector const* b, uint_fast64_t n) const override; - StandardMinMaxLinearEquationSolverSettings const& getSettings() const; - void setSettings(StandardMinMaxLinearEquationSolverSettings const& newSettings); - virtual void clearCache() const override; - virtual ValueType getPrecision() const override; - virtual bool getRelative() const override; - private: - bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - bool solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - bool solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - - bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; - - void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; - - enum class Status { - Converged, TerminatedEarly, MaximalIterationsExceeded, InProgress - }; + protected: // possibly cached data mutable std::unique_ptr> linEqSolverA; mutable std::unique_ptr> auxiliaryRowVector; // A.rowCount() entries - mutable std::unique_ptr> auxiliaryRowGroupVector; // A.rowGroupCount() entries - mutable std::unique_ptr> rowGroupOrdering; // A.rowGroupCount() entries - - Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const; - void reportStatus(Status status, uint64_t iterations) const; - - /// The settings of this solver. - StandardMinMaxLinearEquationSolverSettings settings; /// The factory used to obtain linear equation solvers. std::unique_ptr> linearEquationSolverFactory; @@ -89,44 +40,39 @@ namespace storm { template class StandardMinMaxLinearEquationSolverFactory : public MinMaxLinearEquationSolverFactory { public: - StandardMinMaxLinearEquationSolverFactory(bool trackScheduler = false); - StandardMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, bool trackScheduler = false); - StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, bool trackScheduler = false); + StandardMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); + StandardMinMaxLinearEquationSolverFactory(std::unique_ptr>&& linearEquationSolverFactory, MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); + StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; virtual std::unique_ptr> create(storm::storage::SparseMatrix&& matrix) const override; - StandardMinMaxLinearEquationSolverSettings& getSettings(); - StandardMinMaxLinearEquationSolverSettings const& getSettings() const; - - private: - StandardMinMaxLinearEquationSolverSettings settings; - + protected: std::unique_ptr> linearEquationSolverFactory; }; template class GmmxxMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory { public: - GmmxxMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + GmmxxMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); }; template class EigenMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory { public: - EigenMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + EigenMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); }; template class NativeMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory { public: - NativeMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + NativeMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); }; template class EliminationMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory { public: - EliminationMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + EliminationMinMaxLinearEquationSolverFactory(MinMaxMethodSelection const& method = MinMaxMethodSelection::FROMSETTINGS, bool trackScheduler = false); }; } diff --git a/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp index d47bff510..e320f3497 100644 --- a/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp @@ -458,7 +458,7 @@ namespace storm { } template - TopologicalMinMaxLinearEquationSolverFactory::TopologicalMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { + TopologicalMinMaxLinearEquationSolverFactory::TopologicalMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(MinMaxMethodSelection::Topological, trackScheduler) { // Intentionally left empty. } diff --git a/src/test/storm/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp b/src/test/storm/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp index 8ddc74d27..3d590ef74 100644 --- a/src/test/storm/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp +++ b/src/test/storm/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp @@ -209,8 +209,7 @@ TEST(GmmxxMdpPrctlModelCheckerTest, SchedulerGeneration) { EXPECT_EQ(7ull, mdp->getNumberOfChoices()); - auto solverFactory = std::make_unique>(); - solverFactory->getSettings().setSolutionMethod(storm::solver::StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration); + auto solverFactory = std::make_unique>(storm::solver::MinMaxMethodSelection::PolicyIteration); storm::modelchecker::SparseMdpPrctlModelChecker> checker(*mdp, std::move(solverFactory)); std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("Pmin=? [F \"target\"]"); diff --git a/src/test/storm/solver/GmmxxMinMaxLinearEquationSolverTest.cpp b/src/test/storm/solver/GmmxxMinMaxLinearEquationSolverTest.cpp index b30f329ba..da68456fe 100644 --- a/src/test/storm/solver/GmmxxMinMaxLinearEquationSolverTest.cpp +++ b/src/test/storm/solver/GmmxxMinMaxLinearEquationSolverTest.cpp @@ -109,8 +109,7 @@ TEST(GmmxxMinMaxLinearEquationSolver, SolveWithPolicyIteration) { std::vector x(1); std::vector b = { 0.099, 0.5 }; - auto factory = storm::solver::GmmxxMinMaxLinearEquationSolverFactory(); - factory.getSettings().setSolutionMethod(storm::solver::StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration); + auto factory = storm::solver::GmmxxMinMaxLinearEquationSolverFactory(storm::solver::MinMaxMethodSelection::PolicyIteration); auto solver = factory.create(A); ASSERT_NO_THROW(solver->solveEquations(storm::OptimizationDirection::Minimize, x, b)); diff --git a/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp b/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp index 75f0ddd23..de724e39f 100644 --- a/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp +++ b/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp @@ -80,8 +80,7 @@ TEST(NativeMinMaxLinearEquationSolver, SolveWithPolicyIteration) { std::vector x(1); std::vector b = { 0.099, 0.5 }; - auto factory = storm::solver::NativeMinMaxLinearEquationSolverFactory(); - factory.getSettings().setSolutionMethod(storm::solver::StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration); + auto factory = storm::solver::NativeMinMaxLinearEquationSolverFactory(storm::solver::MinMaxMethodSelection::PolicyIteration); auto solver = factory.create(A); ASSERT_NO_THROW(solver->solveEquations(storm::OptimizationDirection::Minimize, x, b));