diff --git a/src/storm/api/verification.h b/src/storm/api/verification.h index d6d4969d4..868c00be3 100644 --- a/src/storm/api/verification.h +++ b/src/storm/api/verification.h @@ -264,7 +264,7 @@ namespace storm { } template - typename std::enable_if::value, std::unique_ptr>::type verifyWithSparseEngine(storm::Environment const& env, std::shared_ptr> const& mdp, storm::modelchecker::CheckTask const& task) { + typename std::enable_if::value, std::unique_ptr>::type verifyWithSparseEngine(storm::Environment const& env, std::shared_ptr> const& smg, storm::modelchecker::CheckTask const& task) { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Sparse engine cannot verify SMGs with this data type."); } diff --git a/src/storm/logic/FragmentSpecification.cpp b/src/storm/logic/FragmentSpecification.cpp index 96cbaf181..836b020aa 100644 --- a/src/storm/logic/FragmentSpecification.cpp +++ b/src/storm/logic/FragmentSpecification.cpp @@ -63,6 +63,9 @@ namespace storm { rpatl.setLongRunAverageRewardFormulasAllowed(true); rpatl.setLongRunAverageOperatorsAllowed(true); + rpatl.setProbabilityOperatorsAllowed(true); + rpatl.setReachabilityProbabilityFormulasAllowed(true); + return rpatl; } diff --git a/src/storm/logic/PlayerCoalition.cpp b/src/storm/logic/PlayerCoalition.cpp index 316e40de8..45cc683e9 100644 --- a/src/storm/logic/PlayerCoalition.cpp +++ b/src/storm/logic/PlayerCoalition.cpp @@ -8,17 +8,17 @@ namespace storm { PlayerCoalition::PlayerCoalition(std::vector> playerIds) : _playerIds(playerIds) { // Intentionally left empty. } - + std::vector> const& PlayerCoalition::getPlayers() const { return _playerIds; } std::ostream& operator<<(std::ostream& stream, PlayerCoalition const& coalition) { - bool firstItem = true; - for (auto const& id : coalition._playerIds) { - if (firstItem) { firstItem = false; } else { stream << ","; } - stream << id; - } + //bool firstItem = true; + //for (auto const& id : coalition._playerIds) { + // //if(firstItem) { firstItem = false; } else { stream << ","; } + // stream << id; + //} return stream; } } diff --git a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp index e5aeecdd3..239817adc 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp @@ -197,8 +197,6 @@ namespace storm { if (storm::utility::resources::isTerminate()) { break; } - - // If there will be a next iteration, we have to prepare it. if(!gameNondetTs()) { prepareNextIteration(env); diff --git a/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.cpp b/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.cpp index fbfcd796f..c075ab6af 100644 --- a/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.cpp +++ b/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.cpp @@ -10,11 +10,14 @@ #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h" +#include "storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h" +#include "storm/modelchecker/helper/infinitehorizon/SparseNondeterministicGameInfiniteHorizonHelper.h" #include "storm/modelchecker/helper/utility/SetInformationFromCheckTask.h" #include "storm/logic/FragmentSpecification.h" -#include "storm/logic/Coalition.h" +#include "storm/logic/PlayerCoalition.h" #include "storm/storage/BitVector.h" #include "storm/environment/solver/MultiplierEnvironment.h" @@ -23,8 +26,9 @@ #include "storm/settings/modules/GeneralSettings.h" +#include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidPropertyException.h" -#include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/InvalidArgumentException.h" namespace storm { namespace modelchecker { @@ -38,7 +42,6 @@ namespace storm { template bool SparseSmgRpatlModelChecker::canHandleStatic(CheckTask const& checkTask, bool* requiresSingleInitialState) { storm::logic::Formula const& formula = checkTask.getFormula(); - return formula.isInFragment(storm::logic::rpatl()); } @@ -56,17 +59,36 @@ namespace storm { template std::unique_ptr SparseSmgRpatlModelChecker::checkGameFormula(Environment const& env, CheckTask const& checkTask) { + Environment solverEnv = env; + storm::logic::GameFormula const& gameFormula = checkTask.getFormula(); storm::logic::Formula const& subFormula = gameFormula.getSubformula(); + statesOfCoalition = this->getModel().computeStatesOfCoalition(gameFormula.getCoalition()); + if (subFormula.isRewardOperatorFormula()) { return this->checkRewardOperatorFormula(solverEnv, checkTask.substituteFormula(subFormula.asRewardOperatorFormula())); } else if (subFormula.isLongRunAverageOperatorFormula()) { return this->checkLongRunAverageOperatorFormula(solverEnv, checkTask.substituteFormula(subFormula.asLongRunAverageOperatorFormula())); + } else if (subFormula.isProbabilityOperatorFormula()) { + return this->checkProbabilityOperatorFormula(solverEnv, checkTask.substituteFormula(subFormula.asProbabilityOperatorFormula())); } STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Cannot check this property (yet)."); } + template + std::unique_ptr SparseSmgRpatlModelChecker::checkProbabilityOperatorFormula(Environment const& env, CheckTask const& checkTask) { + storm::logic::ProbabilityOperatorFormula const& stateFormula = checkTask.getFormula(); + std::unique_ptr result = this->computeProbabilities(env, checkTask.substituteFormula(stateFormula.getSubformula())); + + if (checkTask.isBoundSet()) { + STORM_LOG_THROW(result->isQuantitative(), storm::exceptions::InvalidOperationException, "Unable to perform comparison operation on non-quantitative result."); + return result->asQuantitativeCheckResult().compareAgainstBound(checkTask.getBoundComparisonType(), checkTask.getBoundThreshold()); + } else { + return result; + } + } + template std::unique_ptr SparseSmgRpatlModelChecker::checkRewardOperatorFormula(Environment const& env, CheckTask const& checkTask) { storm::logic::RewardOperatorFormula const& formula = checkTask.getFormula(); @@ -83,7 +105,14 @@ namespace storm { return result; } - + template + std::unique_ptr SparseSmgRpatlModelChecker::computeProbabilities(Environment const& env, CheckTask const& checkTask) { + storm::logic::Formula const& formula = checkTask.getFormula(); + if (formula.isReachabilityProbabilityFormula()) { + return this->computeReachabilityProbabilities(env, checkTask.substituteFormula(formula.asReachabilityProbabilityFormula())); + } + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << formula << "' is invalid."); + } template std::unique_ptr SparseSmgRpatlModelChecker::computeRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { @@ -94,17 +123,39 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << rewardFormula << "' cannot (yet) be handled."); } + template + std::unique_ptr SparseSmgRpatlModelChecker::computeUntilProbabilities(Environment const& env, CheckTask const& checkTask) { + // Currently we only support computation of reachability probabilities, i.e. the left subformula will always be 'true' (for now). + storm::logic::UntilFormula const& pathFormula = checkTask.getFormula(); + 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::unique_ptr leftResultPointer = this->check(env, pathFormula.getLeftSubformula()); + std::unique_ptr rightResultPointer = this->check(env, pathFormula.getRightSubformula()); + ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); + ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); + storm::solver::SolveGoal foo(this->getModel(), checkTask); + + auto ret = storm::modelchecker::helper::SparseSmgRpatlHelper::computeUntilProbabilities(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), checkTask.isQualitativeSet(), statesOfCoalition, checkTask.isProduceSchedulersSet(), checkTask.getHint()); + std::unique_ptr result(new ExplicitQuantitativeCheckResult(std::move(ret.values))); + if (checkTask.isProduceSchedulersSet() && ret.scheduler) { + result->asExplicitQuantitativeCheckResult().setScheduler(std::move(ret.scheduler)); + } + return result; + } + template std::unique_ptr SparseSmgRpatlModelChecker::computeLongRunAverageProbabilities(Environment const& env, CheckTask const& checkTask) { - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Not implemented."); + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "NYI"); } template std::unique_ptr SparseSmgRpatlModelChecker::computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask); - storm::modelchecker::helper::SparseNondeterministicGameInfiniteHorizonHelper helper(this->getModel().getTransitionMatrix(), this->getModel().getPlayerActionIndices()); + STORM_LOG_THROW(checkTask.isPlayerCoalitionSet(), storm::exceptions::InvalidPropertyException, "No player coalition was set."); + auto coalitionStates = this->getModel().computeStatesOfCoalition(checkTask.getPlayerCoalition()); + std::cout << "Found " << coalitionStates.getNumberOfSetBits() << " states in coalition." << std::endl; + storm::modelchecker::helper::SparseNondeterministicGameInfiniteHorizonHelper helper(this->getModel().getTransitionMatrix(), statesOfCoalition); storm::modelchecker::helper::setInformationFromCheckTaskNondeterministic(helper, checkTask, this->getModel()); - auto values = helper.computeLongRunAverageRewards(env, rewardModel.get()); + auto values = helper.computeLongRunAverageRewards(env, rewardModel.get()); std::unique_ptr result(new ExplicitQuantitativeCheckResult(std::move(values))); if (checkTask.isProduceSchedulersSet()) { diff --git a/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.h b/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.h index abf3243a5..dfbd47762 100644 --- a/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.h +++ b/src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.h @@ -1,11 +1,13 @@ #ifndef STORM_MODELCHECKER_SPARSESMGRPATLMODELCHECKER_H_ #define STORM_MODELCHECKER_SPARSESMGRPATLMODELCHECKER_H_ + #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/models/sparse/Smg.h" #include "storm/utility/solver.h" #include "storm/solver/LinearEquationSolver.h" #include "storm/storage/StronglyConnectedComponent.h" +#include "storm/storage/BitVector.h" namespace storm { namespace modelchecker { @@ -25,11 +27,22 @@ namespace storm { static bool canHandleStatic(CheckTask const& checkTask, bool* requiresSingleInitialState = nullptr); // The implemented methods of the AbstractModelChecker interface. - virtual bool canHandle(CheckTask const& checkTask) const override; - virtual std::unique_ptr checkGameFormula(Environment const& env, CheckTask const& checkTask) override; - - virtual std::unique_ptr computeLongRunAverageProbabilities(Environment const& env, CheckTask const& checkTask) override; - virtual std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; + bool canHandle(CheckTask const& checkTask) const override; + std::unique_ptr checkGameFormula(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr checkProbabilityOperatorFormula(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr checkRewardOperatorFormula(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr checkLongRunAverageOperatorFormula(Environment const& env, CheckTask const& checkTask) override; + + std::unique_ptr computeProbabilities(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr computeRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; + + std::unique_ptr computeUntilProbabilities(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr computeLongRunAverageProbabilities(Environment const& env, CheckTask const& checkTask) override; + std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; + + //void coalitionIndicator(Environment& env, CheckTask const& checkTask); + private: + storm::storage::BitVector statesOfCoalition; }; } // namespace modelchecker } // namespace storm diff --git a/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.cpp b/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.cpp new file mode 100644 index 000000000..9bf5a10e8 --- /dev/null +++ b/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.cpp @@ -0,0 +1,97 @@ +#include "SparseSmgRpatlHelper.h" + +#include "storm/environment/Environment.h" +#include "storm/environment/solver/MultiplierEnvironment.h" +#include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/utility/vector.h" + +#include "storm/modelchecker/rpatl/helper/internal/GameViHelper.h" + + +namespace storm { + namespace modelchecker { + namespace helper { + + template + MDPSparseModelCheckingHelperReturnType SparseSmgRpatlHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::storage::BitVector statesOfCoalition, bool produceScheduler, ModelCheckerHint const& hint) { + // TODO add Kwiatkowska original reference + // TODO FIX solver min max mess + + auto solverEnv = env; + solverEnv.solver().minMax().setMethod(storm::solver::MinMaxMethod::ValueIteration, false); + + // Initialize the solution vector. + std::vector x = std::vector(transitionMatrix.getRowGroupCount() - psiStates.getNumberOfSetBits(), storm::utility::zero()); + std::vector b = transitionMatrix.getConstrainedRowGroupSumVector(~psiStates, psiStates); + + // Reduce matrix to ~Prob1 states- + //STORM_LOG_DEBUG("\n" << transitionMatrix); + storm::storage::SparseMatrix submatrix = transitionMatrix.getSubmatrix(true, ~psiStates, ~psiStates, false); + //STORM_LOG_DEBUG("\n" << submatrix); + + + //STORM_LOG_DEBUG("x = " << storm::utility::vector::toString(x)); + //STORM_LOG_DEBUG("b = " << storm::utility::vector::toString(b)); + + storm::storage::BitVector clippedStatesOfCoalition(statesOfCoalition.size() - psiStates.getNumberOfSetBits()); + //STORM_LOG_DEBUG(psiStates); + //STORM_LOG_DEBUG(statesOfCoalition); + //STORM_LOG_DEBUG(clippedStatesOfCoalition); + + // TODO move this to BitVector-class + auto clippedStatesCounter = 0; + for(uint i = 0; i < psiStates.size(); i++) { + std::cout << i << " : " << psiStates.get(i) << " -> " << statesOfCoalition[i] << std::endl; + if(!psiStates.get(i)) { + clippedStatesOfCoalition.set(clippedStatesCounter, statesOfCoalition[i]); + clippedStatesCounter++; + } + } + //STORM_LOG_DEBUG(clippedStatesOfCoalition); + clippedStatesOfCoalition.complement(); + //STORM_LOG_DEBUG(clippedStatesOfCoalition); + + storm::modelchecker::helper::internal::GameViHelper viHelper(submatrix, clippedStatesOfCoalition); + std::unique_ptr> scheduler; + if (produceScheduler) { + viHelper.setProduceScheduler(true); + } + + viHelper.performValueIteration(env, x, b, goal.direction()); + + STORM_LOG_DEBUG(storm::utility::vector::toString(x)); + if (produceScheduler) { + scheduler = std::make_unique>(expandScheduler(viHelper.extractScheduler(), psiStates)); + STORM_LOG_DEBUG("IsPartial?" << scheduler->isPartialScheduler()); + } + return MDPSparseModelCheckingHelperReturnType(std::move(x), std::move(scheduler)); + } + + template + storm::storage::Scheduler SparseSmgRpatlHelper::expandScheduler(storm::storage::Scheduler scheduler, storm::storage::BitVector psiStates) { + //STORM_LOG_DEBUG(psiStates.size()); + for(uint i = 0; i < 2; i++) { + //STORM_LOG_DEBUG(scheduler.getChoice(i)); + } + storm::storage::Scheduler completeScheduler(psiStates.size()); + uint_fast64_t maybeStatesCounter = 0; + for(uint stateCounter = 0; stateCounter < psiStates.size(); stateCounter++) { + //STORM_LOG_DEBUG(stateCounter << " : " << psiStates.get(stateCounter)); + if(psiStates.get(stateCounter)) { + completeScheduler.setChoice(0, stateCounter); + } else { + completeScheduler.setChoice(scheduler.getChoice(maybeStatesCounter), stateCounter); + maybeStatesCounter++; + } + } + return completeScheduler; + } + + template class SparseSmgRpatlHelper; +#ifdef STORM_HAVE_CARL + template class SparseSmgRpatlHelper; +#endif + } + } +} diff --git a/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h b/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h new file mode 100644 index 000000000..9505bd51f --- /dev/null +++ b/src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h @@ -0,0 +1,47 @@ +#pragma once + +#include + +#include "storm/modelchecker/hints/ModelCheckerHint.h" +#include "storm/modelchecker/prctl/helper/SolutionType.h" +#include "storm/storage/SparseMatrix.h" + +#include "storm/utility/solver.h" +#include "storm/solver/SolveGoal.h" + +#include "storm/modelchecker/prctl/helper/MDPModelCheckingHelperReturnType.h" + +namespace storm { + + class Environment; + + namespace storage { + class BitVector; + } + + namespace models { + namespace sparse { + template + class StandardRewardModel; + } + } + + namespace modelchecker { + class CheckResult; + + namespace helper { + + template + class SparseSmgRpatlHelper { + public: + // TODO should probably be renamed in the future: + + static MDPSparseModelCheckingHelperReturnType computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::storage::BitVector statesOfCoalition, bool produceScheduler, ModelCheckerHint const& hint = ModelCheckerHint()); + + private: + static storm::storage::Scheduler expandScheduler(storm::storage::Scheduler scheduler, storm::storage::BitVector psiStates); + + }; + } + } +} diff --git a/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.cpp b/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.cpp new file mode 100644 index 000000000..4f6fc1e4e --- /dev/null +++ b/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.cpp @@ -0,0 +1,198 @@ +#include "GameViHelper.h" + +#include "storm/environment/Environment.h" +#include "storm/environment/solver/SolverEnvironment.h" +#include "storm/environment/solver/GameSolverEnvironment.h" + + +#include "storm/utility/SignalHandler.h" +#include "storm/utility/vector.h" + +// TODO this will undergo major refactoring as soon as we implement model checking of further properties + +namespace storm { + namespace modelchecker { + namespace helper { + namespace internal { + + template + GameViHelper::GameViHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector statesOfCoalition) : _transitionMatrix(transitionMatrix), _statesOfCoalition(statesOfCoalition) { + } + + template + void GameViHelper::prepareSolversAndMultipliersReachability(const Environment& env) { + // TODO we get whole transitionmatrix and psistates DONE IN smgrpatlhelper + // -> clip statesofcoalition + // -> compute b vector from psiStates + // -> clip transitionmatrix and create multiplier + _multiplier = storm::solver::MultiplierFactory().create(env, _transitionMatrix); + _multiplier->setOptimizationDirectionOverride(_statesOfCoalition); + + _x1IsCurrent = false; + } + + template + void GameViHelper::performValueIteration(Environment const& env, std::vector& x, std::vector b, storm::solver::OptimizationDirection const dir) { + prepareSolversAndMultipliersReachability(env); + ValueType precision = storm::utility::convertNumber(env.solver().game().getPrecision()); + uint64_t maxIter = env.solver().game().getMaximalNumberOfIterations(); + _b = b; + + _x1.assign(_transitionMatrix.getRowGroupCount(), storm::utility::zero()); + _x2 = _x1; + + if (this->isProduceSchedulerSet()) { + if (!this->_producedOptimalChoices.is_initialized()) { + this->_producedOptimalChoices.emplace(); + } + this->_producedOptimalChoices->resize(this->_transitionMatrix.getRowGroupCount()); + } + + uint64_t iter = 0; + std::vector result = x; + while (iter < maxIter) { + ++iter; + performIterationStep(env, dir); + + if (checkConvergence(precision)) { + break; + } + if (storm::utility::resources::isTerminate()) { + break; + } + } + x = xNew(); + + if (isProduceSchedulerSet()) { + // We will be doing one more iteration step and track scheduler choices this time. + performIterationStep(env, dir, &_producedOptimalChoices.get()); + } + } + + template + void GameViHelper::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const dir, std::vector* choices) { + if (!_multiplier) { + prepareSolversAndMultipliersReachability(env); + } + _x1IsCurrent = !_x1IsCurrent; + + // multiplyandreducegaussseidel + // Ax + b + if (choices == nullptr) { + //STORM_LOG_DEBUG("\n" << _transitionMatrix); + //STORM_LOG_DEBUG("xOld = " << storm::utility::vector::toString(xOld())); + //STORM_LOG_DEBUG("b = " << storm::utility::vector::toString(_b)); + //STORM_LOG_DEBUG("xNew = " << storm::utility::vector::toString(xNew())); + _multiplier->multiplyAndReduce(env, dir, xOld(), &_b, xNew()); + //std::cin.get(); + } else { + // Also keep track of the choices made. + _multiplier->multiplyAndReduce(env, dir, xOld(), &_b, xNew(), choices); + } + + // compare applyPointwise + storm::utility::vector::applyPointwise(xOld(), xNew(), xNew(), [&dir] (ValueType const& a, ValueType const& b) -> ValueType { + if(storm::solver::maximize(dir)) { + if(a > b) return a; + else return b; + } else { + if(a > b) return a; + else return b; + } + }); + } + + template + bool GameViHelper::checkConvergence(ValueType threshold) const { + STORM_LOG_ASSERT(_multiplier, "tried to check for convergence without doing an iteration first."); + + // Now check whether the currently produced results are precise enough + STORM_LOG_ASSERT(threshold > storm::utility::zero(), "Did not expect a non-positive threshold."); + auto x1It = xOld().begin(); + auto x1Ite = xOld().end(); + auto x2It = xNew().begin(); + ValueType maxDiff = (*x2It - *x1It); + ValueType minDiff = maxDiff; + // The difference between maxDiff and minDiff is zero at this point. Thus, it doesn't make sense to check the threshold now. + for (++x1It, ++x2It; x1It != x1Ite; ++x1It, ++x2It) { + ValueType diff = (*x2It - *x1It); + // Potentially update maxDiff or minDiff + bool skipCheck = false; + if (maxDiff < diff) { + maxDiff = diff; + } else if (minDiff > diff) { + minDiff = diff; + } else { + skipCheck = true; + } + // Check convergence + if (!skipCheck && (maxDiff - minDiff) > threshold) { + return false; + } + } + // TODO needs checking + return true; + } + + template + void GameViHelper::setProduceScheduler(bool value) { + _produceScheduler = value; + } + + template + bool GameViHelper::isProduceSchedulerSet() const { + return _produceScheduler; + } + + template + std::vector const& GameViHelper::getProducedOptimalChoices() const { + STORM_LOG_ASSERT(this->isProduceSchedulerSet(), "Trying to get the produced optimal choices although no scheduler was requested."); + STORM_LOG_ASSERT(this->_producedOptimalChoices.is_initialized(), "Trying to get the produced optimal choices but none were available. Was there a computation call before?"); + return this->_producedOptimalChoices.get(); + } + + template + std::vector& GameViHelper::getProducedOptimalChoices() { + STORM_LOG_ASSERT(this->isProduceSchedulerSet(), "Trying to get the produced optimal choices although no scheduler was requested."); + STORM_LOG_ASSERT(this->_producedOptimalChoices.is_initialized(), "Trying to get the produced optimal choices but none were available. Was there a computation call before?"); + return this->_producedOptimalChoices.get(); + } + + template + storm::storage::Scheduler GameViHelper::extractScheduler() const{ + auto const& optimalChoices = getProducedOptimalChoices(); + storm::storage::Scheduler scheduler(optimalChoices.size()); + for (uint64_t state = 0; state < optimalChoices.size(); ++state) { + scheduler.setChoice(optimalChoices[state], state); + } + return scheduler; + } + + template + std::vector& GameViHelper::xNew() { + return _x1IsCurrent ? _x1 : _x2; + } + + template + std::vector const& GameViHelper::xNew() const { + return _x1IsCurrent ? _x1 : _x2; + } + + template + std::vector& GameViHelper::xOld() { + return _x1IsCurrent ? _x2 : _x1; + } + + template + std::vector const& GameViHelper::xOld() const { + return _x1IsCurrent ? _x2 : _x1; + } + + template class GameViHelper; +#ifdef STORM_HAVE_CARL + template class GameViHelper; +#endif + } + } + } +} diff --git a/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.h b/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.h new file mode 100644 index 000000000..92f963851 --- /dev/null +++ b/src/storm/modelchecker/rpatl/helper/internal/GameViHelper.h @@ -0,0 +1,77 @@ +#pragma once + +#include "storm/storage/SparseMatrix.h" +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/solver/Multiplier.h" + +namespace storm { + class Environment; + + namespace storage { + template class Scheduler; + } + + namespace modelchecker { + namespace helper { + namespace internal { + + template + class GameViHelper { + public: + GameViHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector statesOfCoalition); + + void prepareSolversAndMultipliersReachability(const Environment& env); + + void performValueIteration(Environment const& env, std::vector& x, std::vector b, storm::solver::OptimizationDirection const dir); + + /*h + * Sets whether an optimal scheduler shall be constructed during the computation + */ + void setProduceScheduler(bool value); + + /*! + * @return whether an optimal scheduler shall be constructed during the computation + */ + bool isProduceSchedulerSet() const; + + storm::storage::Scheduler extractScheduler() const; + private: + void performIterationStep(Environment const& env, storm::solver::OptimizationDirection const dir, std::vector* choices = nullptr); + + /*! + * Checks whether the curently computed value achieves the desired precision + */ + bool checkConvergence(ValueType precision) const; + + std::vector& xNew(); + std::vector const& xNew() const; + + std::vector& xOld(); + std::vector const& xOld() const; + bool _x1IsCurrent; + + /*! + * @pre before calling this, a computation call should have been performed during which scheduler production was enabled. + * @return the produced scheduler of the most recent call. + */ + std::vector const& getProducedOptimalChoices() const; + + /*! + * @pre before calling this, a computation call should have been performed during which scheduler production was enabled. + * @return the produced scheduler of the most recent call. + */ + std::vector& getProducedOptimalChoices(); + + storm::storage::SparseMatrix _transitionMatrix; + storm::storage::BitVector _statesOfCoalition; + std::vector _x1, _x2, _b; + std::unique_ptr> _multiplier; + + bool _produceScheduler = false; + boost::optional> _producedOptimalChoices; + }; + } + } + } +} diff --git a/src/storm/models/sparse/Smg.cpp b/src/storm/models/sparse/Smg.cpp index af9846f19..4e52f40bc 100644 --- a/src/storm/models/sparse/Smg.cpp +++ b/src/storm/models/sparse/Smg.cpp @@ -6,6 +6,7 @@ #include "storm/adapters/RationalFunctionAdapter.h" #include "storm/models/sparse/StandardRewardModel.h" +#include "storm/models/sparse/Mdp.h" #include "storm/exceptions/InvalidArgumentException.h" namespace storm { diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp index a75ab623f..66ae1249a 100644 --- a/src/storm/solver/GmmxxMultiplier.cpp +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -179,7 +179,7 @@ namespace storm { uint64_t selectedChoice; uint64_t currentRow = backwards ? gmmMatrix.nrows() - 1 : 0; - uint64_t currentRowGroup = backwards ? rowGroupIndices.size() - 1 : 0; + uint64_t currentRowGroup = backwards ? rowGroupIndices.size() - 2 : 0; auto row_group_it = backwards ? rowGroupIndices.end() - 2 : rowGroupIndices.begin(); auto row_group_ite = backwards ? rowGroupIndices.begin() - 1 : rowGroupIndices.end() - 1; //if(choices) STORM_LOG_DEBUG(" "); diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 5e71c0503..09bd57726 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -5,6 +5,7 @@ #include "storm/environment/solver/MinMaxSolverEnvironment.h" #include "storm/environment/solver/OviSolverEnvironment.h" +#include "storm/environment/solver/MultiplierEnvironment.h" #include "storm/utility/ConstantsComparator.h" #include "storm/utility/KwekMehlhorn.h" @@ -19,27 +20,27 @@ namespace storm { namespace solver { - + template IterativeMinMaxLinearEquationSolver::IterativeMinMaxLinearEquationSolver(std::unique_ptr>&& linearEquationSolverFactory) : linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { // Intentionally left empty } - + template IterativeMinMaxLinearEquationSolver::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver(A), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { // Intentionally left empty. } - + template IterativeMinMaxLinearEquationSolver::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver(std::move(A)), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { // Intentionally left empty. } - + template MinMaxMethod IterativeMinMaxLinearEquationSolver::getMethod(Environment const& env, bool isExactMode) const { // Adjust the method if none was specified and we want exact or sound computations. auto method = env.solver().minMax().getMethod(); - + if (isExactMode && method != MinMaxMethod::PolicyIteration && method != MinMaxMethod::RationalSearch && method != MinMaxMethod::ViToPi) { if (env.solver().minMax().isMethodSetFromDefault()) { STORM_LOG_INFO("Selecting 'Policy iteration' as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a different method."); @@ -58,7 +59,7 @@ namespace storm { STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::RationalSearch || method == MinMaxMethod::SoundValueIteration || method == MinMaxMethod::IntervalIteration || method == MinMaxMethod::OptimisticValueIteration || method == MinMaxMethod::ViToPi, storm::exceptions::InvalidEnvironmentException, "This solver does not support the selected method."); return method; } - + template bool IterativeMinMaxLinearEquationSolver::internalSolveEquations(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { bool result = false; @@ -87,14 +88,14 @@ namespace storm { default: STORM_LOG_THROW(false, storm::exceptions::InvalidEnvironmentException, "This solver does not implement the selected solution method"); } - + return result; } template bool IterativeMinMaxLinearEquationSolver::solveInducedEquationSystem(Environment const& env, std::unique_ptr>& linearEquationSolver, std::vector const& scheduler, std::vector& x, std::vector& subB, std::vector const& originalB) const { assert(subB.size() == x.size()); - + // Resolve the nondeterminism according to the given scheduler. bool convertToEquationSystem = this->linearEquationSolverFactory->getEquationProblemFormat(env) == LinearEquationSolverProblemFormat::EquationSystem; storm::storage::SparseMatrix submatrix = this->A->selectRowsFromRowGroups(scheduler, convertToEquationSystem); @@ -102,7 +103,7 @@ namespace storm { submatrix.convertToEquationSystem(); } storm::utility::vector::selectVectorValues(subB, scheduler, this->A->getRowGroupIndices(), originalB); - + // Check whether the linear equation solver is already initialized if (!linearEquationSolver) { // Initialize the equation solver @@ -116,14 +117,14 @@ namespace storm { // Solve the equation system for the 'DTMC' and return true upon success return linearEquationSolver->solveEquations(env, x, subB); } - + template bool IterativeMinMaxLinearEquationSolver::solveEquationsPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { // Create the initial scheduler. std::vector scheduler = this->hasInitialScheduler() ? this->getInitialScheduler() : std::vector(this->A->getRowGroupCount()); return performPolicyIteration(env, dir, x, b, std::move(scheduler)); } - + template bool IterativeMinMaxLinearEquationSolver::performPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector&& initialPolicy) const { std::vector scheduler = std::move(initialPolicy); @@ -162,7 +163,7 @@ namespace storm { do { // Solve the equation system for the 'DTMC'. solveInducedEquationSystem(environmentOfSolver, solver, scheduler, x, subB, b); - + // 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) { @@ -172,14 +173,14 @@ namespace storm { 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. @@ -190,12 +191,12 @@ namespace storm { } } } - + // If the scheduler did not improve, we are done. if (!schedulerImproved) { status = SolverStatus::Converged; } - + // Update environment variables. ++iterations; status = this->updateStatus(status, x, dir == storm::OptimizationDirection::Minimize ? SolverGuarantee::GreaterOrEqual : SolverGuarantee::LessOrEqual, iterations, env.solver().minMax().getMaximalNumberOfIterations()); @@ -203,21 +204,21 @@ namespace storm { // Potentially show progress. this->showProgressIterative(iterations); } while (status == SolverStatus::InProgress); - + this->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 == SolverStatus::Converged || status == SolverStatus::TerminatedEarly; } - + template bool IterativeMinMaxLinearEquationSolver::valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const { if (dir == OptimizationDirection::Minimize) { @@ -230,7 +231,7 @@ namespace storm { template MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver::getRequirements(Environment const& env, boost::optional const& direction, bool const& hasInitialScheduler) const { auto method = getMethod(env, storm::NumberTraits::IsExact || env.solver().isForceExact()); - + // Check whether a linear equation solver is needed and potentially start with its requirements bool needsLinEqSolver = false; needsLinEqSolver |= method == MinMaxMethod::PolicyIteration; @@ -259,7 +260,7 @@ namespace storm { requirements.requireUniqueSolution(); } requirements.requireLowerBounds(); - + } else if (method == MinMaxMethod::IntervalIteration) { // Interval iteration requires a unique solution and lower+upper bounds if (!this->hasUniqueSolution()) { @@ -296,18 +297,18 @@ namespace storm { template typename IterativeMinMaxLinearEquationSolver::ValueIterationResult IterativeMinMaxLinearEquationSolver::performValueIteration(Environment const& env, OptimizationDirection dir, std::vector*& currentX, std::vector*& newX, std::vector const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations, uint64_t maximalNumberOfIterations, storm::solver::MultiplicationStyle const& multiplicationStyle) const { - + STORM_LOG_ASSERT(currentX != newX, "Vectors must not be aliased."); - + // Get handle to multiplier. storm::solver::Multiplier const& multiplier = *this->multiplierA; - + // Allow aliased multiplications. bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel; - + // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. uint64_t iterations = currentIterations; - + SolverStatus status = SolverStatus::InProgress; while (status == SolverStatus::InProgress) { // Compute x' = min/max(A*x + b). @@ -318,12 +319,12 @@ namespace storm { } else { multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX); } - + // Determine whether the method converged. if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative)) { status = SolverStatus::Converged; } - + // Update environment variables. std::swap(currentX, newX); ++iterations; @@ -332,7 +333,7 @@ namespace storm { // Potentially show progress. this->showProgressIterative(iterations); } - + return ValueIterationResult(iterations - currentIterations, status); } @@ -367,7 +368,7 @@ namespace storm { } return true; } - + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } @@ -376,14 +377,14 @@ namespace storm { } storm::solver::helper::OptimisticValueIterationHelper helper(*this->A); - + // x has to start with a lower bound. this->createLowerBoundsVector(x); std::vector* lowerX = &x; std::vector* upperX = auxiliaryRowGroupVector.get(); - + auto statusIters = helper.solveEquations(env, lowerX, upperX, b, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision()), @@ -413,14 +414,19 @@ namespace storm { if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } - + + // TODO cleanup + if(env.solver().multiplier().getOptimizationDirectionOverride().is_initialized()) { + multiplierA->setOptimizationDirectionOverride(env.solver().multiplier().getOptimizationDirectionOverride().get()); + } + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - + // By default, we can not provide any guarantee SolverGuarantee guarantee = SolverGuarantee::None; - + if (this->hasInitialScheduler()) { // Solve the equation system induced by the initial scheduler. std::unique_ptr> linEqSolver; @@ -469,7 +475,7 @@ namespace storm { std::vector* newX = auxiliaryRowGroupVector.get(); std::vector* currentX = &x; - + this->startMeasureProgress(); ValueIterationResult result = performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber(env.solver().minMax().getPrecision()), env.solver().minMax().getRelativeTerminationCriterion(), guarantee, 0, env.solver().minMax().getMaximalNumberOfIterations(), env.solver().minMax().getMultiplicationStyle()); @@ -477,27 +483,27 @@ namespace storm { if (currentX == auxiliaryRowGroupVector.get()) { std::swap(x, *currentX); } - + this->reportStatus(result.status, result.iterations); - + // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { this->schedulerChoices = std::vector(this->A->getRowGroupCount()); this->multiplierA->multiplyAndReduce(env, dir, x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get()); } - + if (!this->isCachingEnabled()) { clearCache(); } - + return result.status == SolverStatus::Converged || result.status == SolverStatus::TerminatedEarly; } - + template void preserveOldRelevantValues(std::vector const& allValues, storm::storage::BitVector const& relevantValues, std::vector& oldValues) { storm::utility::vector::selectVectorValues(oldValues, relevantValues, allValues); } - + /*! * This version of value iteration is sound, because it approaches the solution from below and above. This * technique is due to Haddad and Monmege (Interval iteration algorithm for MDPs and IMDPs, TCS 2017) and was @@ -511,28 +517,28 @@ namespace storm { if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } - + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - + // Allow aliased multiplications. bool useGaussSeidelMultiplication = env.solver().minMax().getMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel; - + std::vector* lowerX = &x; this->createLowerBoundsVector(*lowerX); this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount()); std::vector* upperX = this->auxiliaryRowGroupVector.get(); - + std::vector* tmp = nullptr; if (!useGaussSeidelMultiplication) { auxiliaryRowGroupVector2 = std::make_unique>(lowerX->size()); tmp = auxiliaryRowGroupVector2.get(); } - + // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. uint64_t iterations = 0; - + SolverStatus status = SolverStatus::InProgress; bool doConvergenceCheck = true; bool useDiffs = this->hasRelevantValues() && !env.solver().minMax().isSymmetricUpdatesSet(); @@ -636,7 +642,7 @@ namespace storm { status = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, precision, relative) ? SolverStatus::Converged : status; } } - + // Update environment variables. ++iterations; doConvergenceCheck = !doConvergenceCheck; @@ -650,33 +656,33 @@ namespace storm { // Potentially show progress. this->showProgressIterative(iterations); } - + this->reportStatus(status, iterations); // We take the means of the lower and upper bound so we guarantee the desired precision. ValueType two = storm::utility::convertNumber(2.0); storm::utility::vector::applyPointwise(*lowerX, *upperX, *lowerX, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); - + // Since we shuffled the pointer around, we need to write the actual results to the input/output vector x. if (&x == tmp) { std::swap(x, *tmp); } else if (&x == this->auxiliaryRowGroupVector.get()) { std::swap(x, *this->auxiliaryRowGroupVector); } - + // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { this->schedulerChoices = std::vector(this->A->getRowGroupCount()); this->multiplierA->multiplyAndReduce(env, dir, x, &b, *this->auxiliaryRowGroupVector, &this->schedulerChoices.get()); } - + if (!this->isCachingEnabled()) { clearCache(); } - + return status == SolverStatus::Converged; } - + template bool IterativeMinMaxLinearEquationSolver::solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { @@ -690,7 +696,7 @@ namespace storm { } else { this->soundValueIterationHelper = std::make_unique>(std::move(*this->soundValueIterationHelper), x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision())); } - + // Prepare initial bounds for the solution (if given) if (this->hasLowerBound()) { this->soundValueIterationHelper->setLowerBound(this->getLowerBound(true)); @@ -698,16 +704,16 @@ namespace storm { if (this->hasUpperBound()) { this->soundValueIterationHelper->setUpperBound(this->getUpperBound(true)); } - + storm::storage::BitVector const* relevantValuesPtr = nullptr; if (this->hasRelevantValues()) { relevantValuesPtr = &this->getRelevantValues(); } - + SolverStatus status = SolverStatus::InProgress; this->startMeasureProgress(); uint64_t iterations = 0; - + while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) { ++iterations; this->soundValueIterationHelper->performIterationStep(dir, b); @@ -716,12 +722,12 @@ namespace storm { } else { status = this->updateStatus(status, this->hasCustomTerminationCondition() && this->soundValueIterationHelper->checkCustomTerminationCondition(this->getTerminationCondition()), iterations, env.solver().minMax().getMaximalNumberOfIterations()); } - + // Potentially show progress. this->showProgressIterative(iterations); } this->soundValueIterationHelper->setSolutionVector(); - + // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { this->schedulerChoices = std::vector(this->A->getRowGroupCount()); @@ -729,14 +735,14 @@ namespace storm { } this->reportStatus(status, iterations); - + if (!this->isCachingEnabled()) { clearCache(); } - + return status == SolverStatus::Converged; } - + template bool IterativeMinMaxLinearEquationSolver::solveEquationsViToPi(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { // First create an (inprecise) vi solver to get a good initial strategy for the (potentially precise) policy iteration solver. @@ -760,11 +766,11 @@ namespace storm { STORM_LOG_INFO("Found initial policy using Value Iteration. Starting Policy iteration now."); return performPolicyIteration(env, dir, x, b, std::move(initialSched)); } - + template bool IterativeMinMaxLinearEquationSolver::isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix const& matrix, std::vector const& values, std::vector const& b) { storm::utility::ConstantsComparator comparator; - + auto valueIt = values.begin(); auto bIt = b.begin(); for (uint64_t group = 0; group < matrix.getRowGroupCount(); ++group, ++valueIt) { @@ -778,18 +784,18 @@ namespace storm { for (auto endRow = matrix.getRowGroupIndices()[group + 1]; row < endRow; ++row, ++bIt) { ValueType newValue = *bIt; newValue += matrix.multiplyRowWithVector(row, values); - + if ((dir == storm::OptimizationDirection::Minimize && newValue < groupValue) || (dir == storm::OptimizationDirection::Maximize && newValue > groupValue)) { groupValue = newValue; } } - + // If the value does not match the one in the values vector, the given vector is not a solution. if (!comparator.isEqual(groupValue, *valueIt)) { return false; } } - + // Checked all values at this point. return true; } @@ -797,7 +803,7 @@ namespace storm { template template bool IterativeMinMaxLinearEquationSolver::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp) { - + for (uint64_t p = 0; p <= precision; ++p) { storm::utility::kwek_mehlhorn::sharpen(p, x, tmp); @@ -817,18 +823,18 @@ namespace storm { storm::storage::SparseMatrix rationalA = this->A->template toValueType(); std::vector rationalX(x.size()); std::vector rationalB = storm::utility::vector::convertNumericVector(b); - + if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } - + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - + // Forward the call to the core rational search routine. bool converged = solveEquationsRationalSearchHelper(env, dir, *this, rationalA, rationalX, rationalB, *this->A, x, b, *auxiliaryRowGroupVector); - + // Translate back rational result to imprecise result. auto targetIt = x.begin(); for (auto it = rationalX.begin(), ite = rationalX.end(); it != ite; ++it, ++targetIt) { @@ -838,30 +844,30 @@ namespace storm { if (!this->isCachingEnabled()) { this->clearCache(); } - + return converged; } - + template template typename std::enable_if::value && NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { // Version for when the overall value type is exact and the same type is to be used for the imprecise part. - + if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } - + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - + // Forward the call to the core rational search routine. bool converged = solveEquationsRationalSearchHelper(env, dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x); if (!this->isCachingEnabled()) { this->clearCache(); } - + return converged; } @@ -870,10 +876,10 @@ namespace storm { typename std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { // Version for when the overall value type is exact and the imprecise one is not. We first try to solve the // problem using the imprecise data type and fall back to the exact type as needed. - + // Translate A to its imprecise version. storm::storage::SparseMatrix impreciseA = this->A->template toValueType(); - + // Translate x to its imprecise version. std::vector impreciseX(x.size()); { @@ -884,23 +890,23 @@ namespace storm { *targetIt = storm::utility::convertNumber(*sourceIt); } } - + // Create temporary storage for an imprecise x. std::vector impreciseTmpX(x.size()); - + // Translate b to its imprecise version. std::vector impreciseB(b.size()); auto targetIt = impreciseB.begin(); for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) { *targetIt = storm::utility::convertNumber(*sourceIt); } - + // Create imprecise solver from the imprecise data. IterativeMinMaxLinearEquationSolver impreciseSolver(std::make_unique>()); impreciseSolver.setMatrix(impreciseA); impreciseSolver.setCachingEnabled(true); impreciseSolver.multiplierA = storm::solver::MultiplierFactory().create(env, impreciseA); - + bool converged = false; try { // Forward the call to the core rational search routine. @@ -908,7 +914,7 @@ namespace storm { impreciseSolver.clearCache(); } catch (storm::exceptions::PrecisionExceededException const& e) { STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic."); - + if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } @@ -918,7 +924,7 @@ namespace storm { for (auto it = impreciseX.begin(), ite = impreciseX.end(); it != ite; ++it, ++targetIt) { *targetIt = storm::utility::convertNumber(*it); } - + // Get rid of the superfluous data structures. impreciseX = std::vector(); impreciseTmpX = std::vector(); @@ -928,15 +934,15 @@ namespace storm { if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } - + // Forward the call to the core rational search routine, but now with our value type as the imprecise value type. converged = solveEquationsRationalSearchHelper(env, dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x); } - + if (!this->isCachingEnabled()) { this->clearCache(); } - + return converged; } @@ -945,12 +951,12 @@ namespace storm { static std::vector* getTemporary(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { return &rationalX; } - + static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { // Nothing to do. } }; - + template struct TemporaryHelper { static std::vector* getTemporary(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { @@ -960,7 +966,7 @@ namespace storm { static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { if (&rationalX == rationalSolution) { // In this case, the rational solution is in place. - + // However, since the rational solution is no alias to current x, the imprecise solution is stored // in current x and and rational x is not an alias to x, we can swap the contents of currentX to x. std::swap(x, *currentX); @@ -971,7 +977,7 @@ namespace storm { } } }; - + template template bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, IterativeMinMaxLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const { @@ -989,29 +995,29 @@ namespace storm { while (status == SolverStatus::InProgress && overallIterations < env.solver().minMax().getMaximalNumberOfIterations()) { // Perform value iteration with the current precision. typename IterativeMinMaxLinearEquationSolver::ValueIterationResult result = impreciseSolver.performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber(precision), env.solver().minMax().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations, env.solver().minMax().getMaximalNumberOfIterations(), env.solver().minMax().getMultiplicationStyle()); - + // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x. - + ++valueIterationInvocations; STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations."); - + // Count the iterations. overallIterations += result.iterations; - + // Compute maximal precision until which to sharpen. uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); - + // Make sure that currentX and rationalX are not aliased. std::vector* temporaryRational = TemporaryHelper::getTemporary(rationalX, currentX, newX); - + // Sharpen solution and place it in the temporary rational. bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, *temporaryRational); - + // After sharpen, if a solution was found, it is contained in the free rational. - + if (foundSolution) { status = SolverStatus::Converged; - + TemporaryHelper::swapSolutions(rationalX, temporaryRational, x, currentX, newX); } else { // Increase the precision. @@ -1020,14 +1026,14 @@ namespace storm { status = this->updateStatus(status, false, overallIterations, env.solver().minMax().getMaximalNumberOfIterations()); } - + // Swap the two vectors if the current result is not in the original x. if (currentX != originalX) { std::swap(x, tmpX); } - + this->reportStatus(status, overallIterations); - + return status == SolverStatus::Converged || status == SolverStatus::TerminatedEarly; } @@ -1035,18 +1041,18 @@ namespace storm { bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearch(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { return solveEquationsRationalSearchHelper(env, dir, x, b); } - + 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)) { @@ -1065,7 +1071,7 @@ namespace storm { *choice = optimalRow - this->A->getRowGroupIndices()[group]; } } - + template void IterativeMinMaxLinearEquationSolver::clearCache() const { multiplierA.reset(); @@ -1075,9 +1081,9 @@ namespace storm { optimisticValueIterationHelper.reset(); StandardMinMaxLinearEquationSolver::clearCache(); } - + template class IterativeMinMaxLinearEquationSolver; - + #ifdef STORM_HAVE_CARL template class IterativeMinMaxLinearEquationSolver; #endif diff --git a/src/storm/utility/Engine.cpp b/src/storm/utility/Engine.cpp index e98ccfea2..adb9ca6e9 100644 --- a/src/storm/utility/Engine.cpp +++ b/src/storm/utility/Engine.cpp @@ -26,7 +26,7 @@ namespace storm { namespace utility { - + // Returns a list of all available engines (excluding Unknown) std::vector getEngines() { std::vector res; @@ -35,7 +35,7 @@ namespace storm { } return res; } - + std::string toString(Engine const& engine) { switch (engine) { case Engine::Sparse: @@ -61,7 +61,7 @@ namespace storm { return "UNKNOWN"; } } - + std::ostream& operator<<(std::ostream& os, Engine const& engine) { os << toString(engine); return os; @@ -80,7 +80,7 @@ namespace storm { STORM_LOG_ERROR("The engine '" << engineStr << "' was not found."); return Engine::Unknown; } - + storm::builder::BuilderType getBuilderType(Engine const& engine) { switch (engine) { case Engine::Sparse: @@ -217,7 +217,7 @@ namespace storm { STORM_LOG_ERROR("The selected combination of engine (" << engine << ") and model type (" << modelType << ") does not seem to be supported for this value type."); return false; } - + template bool canHandle(storm::utility::Engine const& engine, std::vector const& properties, storm::storage::SymbolicModelDescription const& modelDescription) { // Check handability of properties based on model type @@ -233,7 +233,7 @@ namespace storm { // Check whether the model builder can handle the model description return storm::builder::canHandle(getBuilderType(engine), modelDescription, properties); } - + // explicit template instantiations. template bool canHandle(storm::utility::Engine const&, std::vector const&, storm::storage::SymbolicModelDescription const&); template bool canHandle(storm::utility::Engine const&, std::vector const&, storm::storage::SymbolicModelDescription const&);