Browse Source

Merge pull request 'Merge simple reachability for SMG' (#5) from smg_reachability into main

Reviewed-on: http://git.pranger.xyz/TEMPEST/tempest-devel/pulls/5

Going to receive a major cleanup just as #4
tempestpy_adaptions
Stefan Pranger 4 years ago
parent
commit
f3a2f89b7a
  1. 2
      src/storm/api/verification.h
  2. 3
      src/storm/logic/FragmentSpecification.cpp
  3. 12
      src/storm/logic/PlayerCoalition.cpp
  4. 2
      src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp
  5. 65
      src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.cpp
  6. 23
      src/storm/modelchecker/rpatl/SparseSmgRpatlModelChecker.h
  7. 97
      src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.cpp
  8. 47
      src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h
  9. 198
      src/storm/modelchecker/rpatl/helper/internal/GameViHelper.cpp
  10. 77
      src/storm/modelchecker/rpatl/helper/internal/GameViHelper.h
  11. 1
      src/storm/models/sparse/Smg.cpp
  12. 2
      src/storm/solver/GmmxxMultiplier.cpp
  13. 234
      src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
  14. 12
      src/storm/utility/Engine.cpp

2
src/storm/api/verification.h

@ -264,7 +264,7 @@ namespace storm {
}
template<typename ValueType>
typename std::enable_if<std::is_same<ValueType, storm::RationalFunction>::value, std::unique_ptr<storm::modelchecker::CheckResult>>::type verifyWithSparseEngine(storm::Environment const& env, std::shared_ptr<storm::models::sparse::Smg<ValueType>> const& mdp, storm::modelchecker::CheckTask<storm::logic::Formula, ValueType> const& task) {
typename std::enable_if<std::is_same<ValueType, storm::RationalFunction>::value, std::unique_ptr<storm::modelchecker::CheckResult>>::type verifyWithSparseEngine(storm::Environment const& env, std::shared_ptr<storm::models::sparse::Smg<ValueType>> const& smg, storm::modelchecker::CheckTask<storm::logic::Formula, ValueType> const& task) {
STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Sparse engine cannot verify SMGs with this data type.");
}

3
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;
}

12
src/storm/logic/PlayerCoalition.cpp

@ -8,17 +8,17 @@ namespace storm {
PlayerCoalition::PlayerCoalition(std::vector<boost::variant<std::string, storm::storage::PlayerIndex>> playerIds) : _playerIds(playerIds) {
// Intentionally left empty.
}
std::vector<boost::variant<std::string, storm::storage::PlayerIndex>> 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;
}
}

2
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);

65
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<typename SparseSmgModelType>
bool SparseSmgRpatlModelChecker<SparseSmgModelType>::canHandleStatic(CheckTask<storm::logic::Formula, ValueType> const& checkTask, bool* requiresSingleInitialState) {
storm::logic::Formula const& formula = checkTask.getFormula();
return formula.isInFragment(storm::logic::rpatl());
}
@ -56,17 +59,36 @@ namespace storm {
template<typename SparseSmgModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<SparseSmgModelType>::checkGameFormula(Environment const& env, CheckTask<storm::logic::GameFormula, ValueType> 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<typename ModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<ModelType>::checkProbabilityOperatorFormula(Environment const& env, CheckTask<storm::logic::ProbabilityOperatorFormula, ValueType> const& checkTask) {
storm::logic::ProbabilityOperatorFormula const& stateFormula = checkTask.getFormula();
std::unique_ptr<CheckResult> 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<ValueType>().compareAgainstBound(checkTask.getBoundComparisonType(), checkTask.getBoundThreshold());
} else {
return result;
}
}
template<typename ModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<ModelType>::checkRewardOperatorFormula(Environment const& env, CheckTask<storm::logic::RewardOperatorFormula, ValueType> const& checkTask) {
storm::logic::RewardOperatorFormula const& formula = checkTask.getFormula();
@ -83,7 +105,14 @@ namespace storm {
return result;
}
template<typename ModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<ModelType>::computeProbabilities(Environment const& env, CheckTask<storm::logic::Formula, ValueType> 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<typename ModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<ModelType>::computeRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::Formula, ValueType> const& checkTask) {
@ -94,17 +123,39 @@ namespace storm {
STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << rewardFormula << "' cannot (yet) be handled.");
}
template<typename ModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<ModelType>::computeUntilProbabilities(Environment const& env, CheckTask<storm::logic::UntilFormula, ValueType> 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<CheckResult> leftResultPointer = this->check(env, pathFormula.getLeftSubformula());
std::unique_ptr<CheckResult> rightResultPointer = this->check(env, pathFormula.getRightSubformula());
ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
storm::solver::SolveGoal<ValueType> foo(this->getModel(), checkTask);
auto ret = storm::modelchecker::helper::SparseSmgRpatlHelper<ValueType>::computeUntilProbabilities(env, storm::solver::SolveGoal<ValueType>(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), checkTask.isQualitativeSet(), statesOfCoalition, checkTask.isProduceSchedulersSet(), checkTask.getHint());
std::unique_ptr<CheckResult> result(new ExplicitQuantitativeCheckResult<ValueType>(std::move(ret.values)));
if (checkTask.isProduceSchedulersSet() && ret.scheduler) {
result->asExplicitQuantitativeCheckResult<ValueType>().setScheduler(std::move(ret.scheduler));
}
return result;
}
template<typename SparseSmgModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<SparseSmgModelType>::computeLongRunAverageProbabilities(Environment const& env, CheckTask<storm::logic::StateFormula, ValueType> const& checkTask) {
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Not implemented.");
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "NYI");
}
template<typename SparseSmgModelType>
std::unique_ptr<CheckResult> SparseSmgRpatlModelChecker<SparseSmgModelType>::computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) {
auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask);
storm::modelchecker::helper::SparseNondeterministicGameInfiniteHorizonHelper<ValueType> 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<ValueType> 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<CheckResult> result(new ExplicitQuantitativeCheckResult<ValueType>(std::move(values)));
if (checkTask.isProduceSchedulersSet()) {

23
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<storm::logic::Formula, ValueType> const& checkTask, bool* requiresSingleInitialState = nullptr);
// The implemented methods of the AbstractModelChecker interface.
virtual bool canHandle(CheckTask<storm::logic::Formula, ValueType> const& checkTask) const override;
virtual std::unique_ptr<CheckResult> checkGameFormula(Environment const& env, CheckTask<storm::logic::GameFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeLongRunAverageProbabilities(Environment const& env, CheckTask<storm::logic::StateFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
bool canHandle(CheckTask<storm::logic::Formula, ValueType> const& checkTask) const override;
std::unique_ptr<CheckResult> checkGameFormula(Environment const& env, CheckTask<storm::logic::GameFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> checkProbabilityOperatorFormula(Environment const& env, CheckTask<storm::logic::ProbabilityOperatorFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> checkRewardOperatorFormula(Environment const& env, CheckTask<storm::logic::RewardOperatorFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> checkLongRunAverageOperatorFormula(Environment const& env, CheckTask<storm::logic::LongRunAverageOperatorFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> computeProbabilities(Environment const& env, CheckTask<storm::logic::Formula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> computeRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::Formula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> computeUntilProbabilities(Environment const& env, CheckTask<storm::logic::UntilFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> computeLongRunAverageProbabilities(Environment const& env, CheckTask<storm::logic::StateFormula, ValueType> const& checkTask) override;
std::unique_ptr<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
//void coalitionIndicator(Environment& env, CheckTask<storm::logic::GameFormula, ValueType> const& checkTask);
private:
storm::storage::BitVector statesOfCoalition;
};
} // namespace modelchecker
} // namespace storm

97
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<typename ValueType>
MDPSparseModelCheckingHelperReturnType<ValueType> SparseSmgRpatlHelper<ValueType>::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> 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<ValueType> x = std::vector<ValueType>(transitionMatrix.getRowGroupCount() - psiStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
std::vector<ValueType> b = transitionMatrix.getConstrainedRowGroupSumVector(~psiStates, psiStates);
// Reduce matrix to ~Prob1 states-
//STORM_LOG_DEBUG("\n" << transitionMatrix);
storm::storage::SparseMatrix<ValueType> 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<ValueType> viHelper(submatrix, clippedStatesOfCoalition);
std::unique_ptr<storm::storage::Scheduler<ValueType>> 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<storm::storage::Scheduler<ValueType>>(expandScheduler(viHelper.extractScheduler(), psiStates));
STORM_LOG_DEBUG("IsPartial?" << scheduler->isPartialScheduler());
}
return MDPSparseModelCheckingHelperReturnType<ValueType>(std::move(x), std::move(scheduler));
}
template<typename ValueType>
storm::storage::Scheduler<ValueType> SparseSmgRpatlHelper<ValueType>::expandScheduler(storm::storage::Scheduler<ValueType> 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<ValueType> 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<double>;
#ifdef STORM_HAVE_CARL
template class SparseSmgRpatlHelper<storm::RationalNumber>;
#endif
}
}
}

47
src/storm/modelchecker/rpatl/helper/SparseSmgRpatlHelper.h

@ -0,0 +1,47 @@
#pragma once
#include <vector>
#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 <typename ValueType>
class StandardRewardModel;
}
}
namespace modelchecker {
class CheckResult;
namespace helper {
template <typename ValueType>
class SparseSmgRpatlHelper {
public:
// TODO should probably be renamed in the future:
static MDPSparseModelCheckingHelperReturnType<ValueType> computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> 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<ValueType> expandScheduler(storm::storage::Scheduler<ValueType> scheduler, storm::storage::BitVector psiStates);
};
}
}
}

198
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 <typename ValueType>
GameViHelper<ValueType>::GameViHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector statesOfCoalition) : _transitionMatrix(transitionMatrix), _statesOfCoalition(statesOfCoalition) {
}
template <typename ValueType>
void GameViHelper<ValueType>::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<ValueType>().create(env, _transitionMatrix);
_multiplier->setOptimizationDirectionOverride(_statesOfCoalition);
_x1IsCurrent = false;
}
template <typename ValueType>
void GameViHelper<ValueType>::performValueIteration(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> b, storm::solver::OptimizationDirection const dir) {
prepareSolversAndMultipliersReachability(env);
ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().game().getPrecision());
uint64_t maxIter = env.solver().game().getMaximalNumberOfIterations();
_b = b;
_x1.assign(_transitionMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
_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<ValueType> 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 <typename ValueType>
void GameViHelper<ValueType>::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const dir, std::vector<uint64_t>* 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<ValueType, ValueType, ValueType>(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 <typename ValueType>
bool GameViHelper<ValueType>::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<ValueType>(), "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 <typename ValueType>
void GameViHelper<ValueType>::setProduceScheduler(bool value) {
_produceScheduler = value;
}
template <typename ValueType>
bool GameViHelper<ValueType>::isProduceSchedulerSet() const {
return _produceScheduler;
}
template <typename ValueType>
std::vector<uint64_t> const& GameViHelper<ValueType>::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 <typename ValueType>
std::vector<uint64_t>& GameViHelper<ValueType>::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 <typename ValueType>
storm::storage::Scheduler<ValueType> GameViHelper<ValueType>::extractScheduler() const{
auto const& optimalChoices = getProducedOptimalChoices();
storm::storage::Scheduler<ValueType> scheduler(optimalChoices.size());
for (uint64_t state = 0; state < optimalChoices.size(); ++state) {
scheduler.setChoice(optimalChoices[state], state);
}
return scheduler;
}
template <typename ValueType>
std::vector<ValueType>& GameViHelper<ValueType>::xNew() {
return _x1IsCurrent ? _x1 : _x2;
}
template <typename ValueType>
std::vector<ValueType> const& GameViHelper<ValueType>::xNew() const {
return _x1IsCurrent ? _x1 : _x2;
}
template <typename ValueType>
std::vector<ValueType>& GameViHelper<ValueType>::xOld() {
return _x1IsCurrent ? _x2 : _x1;
}
template <typename ValueType>
std::vector<ValueType> const& GameViHelper<ValueType>::xOld() const {
return _x1IsCurrent ? _x2 : _x1;
}
template class GameViHelper<double>;
#ifdef STORM_HAVE_CARL
template class GameViHelper<storm::RationalNumber>;
#endif
}
}
}
}

77
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 <typename VT> class Scheduler;
}
namespace modelchecker {
namespace helper {
namespace internal {
template <typename ValueType>
class GameViHelper {
public:
GameViHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector statesOfCoalition);
void prepareSolversAndMultipliersReachability(const Environment& env);
void performValueIteration(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> 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<ValueType> extractScheduler() const;
private:
void performIterationStep(Environment const& env, storm::solver::OptimizationDirection const dir, std::vector<uint64_t>* choices = nullptr);
/*!
* Checks whether the curently computed value achieves the desired precision
*/
bool checkConvergence(ValueType precision) const;
std::vector<ValueType>& xNew();
std::vector<ValueType> const& xNew() const;
std::vector<ValueType>& xOld();
std::vector<ValueType> 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<uint64_t> 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<uint64_t>& getProducedOptimalChoices();
storm::storage::SparseMatrix<ValueType> _transitionMatrix;
storm::storage::BitVector _statesOfCoalition;
std::vector<ValueType> _x1, _x2, _b;
std::unique_ptr<storm::solver::Multiplier<ValueType>> _multiplier;
bool _produceScheduler = false;
boost::optional<std::vector<uint64_t>> _producedOptimalChoices;
};
}
}
}
}

1
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 {

2
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(" ");

234
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<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty
}
template<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(A), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty.
}
template<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(std::move(A)), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty.
}
template<typename ValueType>
MinMaxMethod IterativeMinMaxLinearEquationSolver<ValueType>::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<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::internalSolveEquations(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> 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<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveInducedEquationSystem(Environment const& env, std::unique_ptr<LinearEquationSolver<ValueType>>& linearEquationSolver, std::vector<uint64_t> const& scheduler, std::vector<ValueType>& x, std::vector<ValueType>& subB, std::vector<ValueType> 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<ValueType> submatrix = this->A->selectRowsFromRowGroups(scheduler, convertToEquationSystem);
@ -102,7 +103,7 @@ namespace storm {
submatrix.convertToEquationSystem();
}
storm::utility::vector::selectVectorValues<ValueType>(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<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Create the initial scheduler.
std::vector<storm::storage::sparse::state_type> scheduler = this->hasInitialScheduler() ? this->getInitialScheduler() : std::vector<storm::storage::sparse::state_type>(this->A->getRowGroupCount());
return performPolicyIteration(env, dir, x, b, std::move(scheduler));
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::performPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<storm::storage::sparse::state_type>&& initialPolicy) const {
std::vector<storm::storage::sparse::state_type> 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<ValueType>();
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<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const {
if (dir == OptimizationDirection::Minimize) {
@ -230,7 +231,7 @@ namespace storm {
template<typename ValueType>
MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver<ValueType>::getRequirements(Environment const& env, boost::optional<storm::solver::OptimizationDirection> const& direction, bool const& hasInitialScheduler) const {
auto method = getMethod(env, storm::NumberTraits<ValueType>::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 ValueType>
typename IterativeMinMaxLinearEquationSolver<ValueType>::ValueIterationResult IterativeMinMaxLinearEquationSolver<ValueType>::performValueIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> 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<ValueType> 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<ValueType>(*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<std::vector<ValueType>>(this->A->getRowGroupCount());
}
@ -376,14 +377,14 @@ namespace storm {
}
storm::solver::helper::OptimisticValueIterationHelper<ValueType> helper(*this->A);
// x has to start with a lower bound.
this->createLowerBoundsVector(x);
std::vector<ValueType>* lowerX = &x;
std::vector<ValueType>* upperX = auxiliaryRowGroupVector.get();
auto statusIters = helper.solveEquations(env, lowerX, upperX, b,
env.solver().minMax().getRelativeTerminationCriterion(),
storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()),
@ -413,14 +414,19 @@ namespace storm {
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().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<std::vector<ValueType>>(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<storm::solver::LinearEquationSolver<ValueType>> linEqSolver;
@ -469,7 +475,7 @@ namespace storm {
std::vector<ValueType>* newX = auxiliaryRowGroupVector.get();
std::vector<ValueType>* currentX = &x;
this->startMeasureProgress();
ValueIterationResult result = performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber<ValueType>(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<uint_fast64_t>(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<typename ValueType>
void preserveOldRelevantValues(std::vector<ValueType> const& allValues, storm::storage::BitVector const& relevantValues, std::vector<ValueType>& 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<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Allow aliased multiplications.
bool useGaussSeidelMultiplication = env.solver().minMax().getMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;
std::vector<ValueType>* lowerX = &x;
this->createLowerBoundsVector(*lowerX);
this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount());
std::vector<ValueType>* upperX = this->auxiliaryRowGroupVector.get();
std::vector<ValueType>* tmp = nullptr;
if (!useGaussSeidelMultiplication) {
auxiliaryRowGroupVector2 = std::make_unique<std::vector<ValueType>>(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<ValueType>(*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<ValueType>(2.0);
storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(*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<uint_fast64_t>(this->A->getRowGroupCount());
this->multiplierA->multiplyAndReduce(env, dir, x, &b, *this->auxiliaryRowGroupVector, &this->schedulerChoices.get());
}
if (!this->isCachingEnabled()) {
clearCache();
}
return status == SolverStatus::Converged;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
@ -690,7 +696,7 @@ namespace storm {
} else {
this->soundValueIterationHelper = std::make_unique<storm::solver::helper::SoundValueIterationHelper<ValueType>>(std::move(*this->soundValueIterationHelper), x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(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<uint_fast64_t>(this->A->getRowGroupCount());
@ -729,14 +735,14 @@ namespace storm {
}
this->reportStatus(status, iterations);
if (!this->isCachingEnabled()) {
clearCache();
}
return status == SolverStatus::Converged;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsViToPi(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> 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<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& values, std::vector<ValueType> const& b) {
storm::utility::ConstantsComparator<ValueType> 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<typename ValueType>
template<typename RationalType, typename ImpreciseType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& 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<storm::RationalNumber> rationalA = this->A->template toValueType<storm::RationalNumber>();
std::vector<storm::RationalNumber> rationalX(x.size());
std::vector<storm::RationalNumber> rationalB = storm::utility::vector::convertNumericVector<storm::RationalNumber>(b);
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Forward the call to the core rational search routine.
bool converged = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(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<typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> 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<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Forward the call to the core rational search routine.
bool converged = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(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<!std::is_same<ValueType, ImpreciseType>::value, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> 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<ImpreciseType> impreciseA = this->A->template toValueType<ImpreciseType>();
// Translate x to its imprecise version.
std::vector<ImpreciseType> impreciseX(x.size());
{
@ -884,23 +890,23 @@ namespace storm {
*targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
}
// Create temporary storage for an imprecise x.
std::vector<ImpreciseType> impreciseTmpX(x.size());
// Translate b to its imprecise version.
std::vector<ImpreciseType> impreciseB(b.size());
auto targetIt = impreciseB.begin();
for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) {
*targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
// Create imprecise solver from the imprecise data.
IterativeMinMaxLinearEquationSolver<ImpreciseType> impreciseSolver(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ImpreciseType>>());
impreciseSolver.setMatrix(impreciseA);
impreciseSolver.setCachingEnabled(true);
impreciseSolver.multiplierA = storm::solver::MultiplierFactory<ImpreciseType>().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<std::vector<ValueType>>(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<ValueType>(*it);
}
// Get rid of the superfluous data structures.
impreciseX = std::vector<ImpreciseType>();
impreciseTmpX = std::vector<ImpreciseType>();
@ -928,15 +934,15 @@ namespace storm {
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().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<ValueType, ValueType>(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<RationalType>* getTemporary(std::vector<RationalType>& rationalX, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
return &rationalX;
}
static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<ImpreciseType>& x, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
// Nothing to do.
}
};
template<typename RationalType>
struct TemporaryHelper<RationalType, RationalType> {
static std::vector<RationalType>* getTemporary(std::vector<RationalType>& rationalX, std::vector<RationalType>*& currentX, std::vector<RationalType>*& newX) {
@ -960,7 +966,7 @@ namespace storm {
static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<RationalType>& x, std::vector<RationalType>*& currentX, std::vector<RationalType>*& 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<typename ValueType>
template<typename RationalType, typename ImpreciseType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, IterativeMinMaxLinearEquationSolver<ImpreciseType> const& impreciseSolver, storm::storage::SparseMatrix<RationalType> const& rationalA, std::vector<RationalType>& rationalX, std::vector<RationalType> const& rationalB, storm::storage::SparseMatrix<ImpreciseType> const& A, std::vector<ImpreciseType>& x, std::vector<ImpreciseType> const& b, std::vector<ImpreciseType>& 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<ImpreciseType>::ValueIterationResult result = impreciseSolver.performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(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<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
// Make sure that currentX and rationalX are not aliased.
std::vector<RationalType>* temporaryRational = TemporaryHelper<RationalType, ImpreciseType>::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<RationalType, ImpreciseType>::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<ValueType>::solveEquationsRationalSearch(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
return solveEquationsRationalSearchHelper<double>(env, dir, x, b);
}
template<typename ValueType>
void IterativeMinMaxLinearEquationSolver<ValueType>::computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> 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<typename ValueType>
void IterativeMinMaxLinearEquationSolver<ValueType>::clearCache() const {
multiplierA.reset();
@ -1075,9 +1081,9 @@ namespace storm {
optimisticValueIterationHelper.reset();
StandardMinMaxLinearEquationSolver<ValueType>::clearCache();
}
template class IterativeMinMaxLinearEquationSolver<double>;
#ifdef STORM_HAVE_CARL
template class IterativeMinMaxLinearEquationSolver<storm::RationalNumber>;
#endif

12
src/storm/utility/Engine.cpp

@ -26,7 +26,7 @@
namespace storm {
namespace utility {
// Returns a list of all available engines (excluding Unknown)
std::vector<Engine> getEngines() {
std::vector<Engine> 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 <typename ValueType>
bool canHandle(storm::utility::Engine const& engine, std::vector<storm::jani::Property> 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<ValueType>(getBuilderType(engine), modelDescription, properties);
}
// explicit template instantiations.
template bool canHandle<double>(storm::utility::Engine const&, std::vector<storm::jani::Property> const&, storm::storage::SymbolicModelDescription const&);
template bool canHandle<storm::RationalNumber>(storm::utility::Engine const&, std::vector<storm::jani::Property> const&, storm::storage::SymbolicModelDescription const&);

Loading…
Cancel
Save