diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp new file mode 100644 index 000000000..d6dc8eb3c --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp @@ -0,0 +1,209 @@ +#include "SparseDeterministicInfiniteHorizonHelper.h" + +#include "storm/modelchecker/helper/infinitehorizon/internal/ComponentUtility.h" +#include "storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h" + + +#include "storm/storage/SparseMatrix.h" +#include "storm/storage/StronglyConnectedComponentDecomposition.h" +#include "storm/storage/Scheduler.h" + +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/Multiplier.h" +#include "storm/solver/LpSolver.h" + +#include "storm/utility/SignalHandler.h" +#include "storm/utility/solver.h" +#include "storm/utility/vector.h" + +#include "storm/environment/solver/LongRunAverageSolverEnvironment.h" +#include "storm/environment/solver/MinMaxSolverEnvironment.h" + +#include "storm/exceptions/UnmetRequirementException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + SparseDeterministicInfiniteHorizonHelper::SparseDeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix) : SparseInfiniteHorizonHelper(transitionMatrix) { + // Intentionally left empty. + } + + template + SparseDeterministicInfiniteHorizonHelper::SparseDeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates) : SparseInfiniteHorizonHelper(transitionMatrix, exitRates) { + // For the CTMC case we assert that the caller actually provided the probabilistic transitions + STORM_LOG_ASSERT(this->_transitionMatrix.isProbabilistic(), "Non-probabilistic transitions"); + } + + template + void SparseDeterministicInfiniteHorizonHelper::createDecomposition() { + if (this->_longRunComponentDecomposition == nullptr) { + // The decomposition has not been provided or computed, yet. + this->_computedLongRunComponentDecomposition = std::make_unique>(this->_transitionMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs()); + this->_longRunComponentDecomposition = this->_computedLongRunComponentDecomposition.get(); + } + } + + template + ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) { + // For deterministic models, we compute the LRA for a BSCC + + STORM_LOG_ASSERT(!this->isProduceSchedulerSet(), "Scheduler production enabled for deterministic model."); + + auto trivialResult = computeLraForTrivialBscc(env, stateRewardsGetter, actionRewardsGetter, component); + if (trivialResult.first) { + return trivialResult.second; + } + + // Solve nontrivial BSCC with the method specified in the settings + storm::solver::LraMethod method = env.solver().lra().getDetLraMethod(); + if ((storm::NumberTraits::IsExact || env.solver().isForceExact()) && env.solver().lra().isDetLraMethodSetFromDefault() && method == storm::solver::LraMethod::ValueIteration) { + method = storm::solver::LraMethod::GainBiasEquations; + STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method."); + } else if (env.solver().isForceSoundness() && env.solver().lra().isDetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { + method = storm::solver::LraMethod::ValueIteration; + STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method."); + } + STORM_LOG_TRACE("Computing LRA for BSCC of size " << component.size() << " using '" << storm::solver::toString(method) << "'."); + if (method == storm::solver::LraMethod::ValueIteration) { + return computeLraForBsccVi(env, stateRewardsGetter, actionRewardsGetter, component); + }/* else if (method == storm::solver::LraMethod::LraDistributionEquations) { + // We only need the first element of the pair as the lra distribution is not relevant at this point. + return computeLongRunAveragesForBsccLraDistr(env, bscc, rateMatrix, valueGetter, exitRateVector).first; + } + STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to " << storm::solver::toString(storm::solver::LraMethod::GainBiasEquations) << "."); + // We don't need the bias values + return computeLongRunAveragesForBsccGainBias(env, bscc, rateMatrix, valueGetter, exitRateVector).first;*/ + } + + template + std::pair SparseDeterministicInfiniteHorizonHelper::computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) { + + // For deterministic models, we can catch the case where all values are the same. This includes the special case where the BSCC consist only of just one state. + bool first = true; + ValueType val = storm::utility::zero(); + for (auto const& element : component) { + auto state = internal::getComponentElementState(element); + STORM_LOG_ASSERT(state == *internal::getComponentElementChoicesBegin(element), "Unexpected choice index at state " << state << " of deterministic model."); + ValueType curr = stateRewardsGetter(state) + (this->isContinuousTime() ? (*this->_exitRates)[state] * actionRewardsGetter(state) : actionRewardsGetter(state)); + if (first) { + first = false; + } else if (val != curr) { + return {false, storm::utility::zero()}; + } + } + // All values are the same + return {true, val}; + } + + + template + ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForBsccVi(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& bscc) { + + // Collect parameters of the computation + ValueType aperiodicFactor = storm::utility::convertNumber(env.solver().lra().getAperiodicFactor()); + + // Now create a helper and perform the algorithm + if (this->isContinuousTime()) { + // We assume a CTMC (with deterministic timed states and no instant states) + storm::modelchecker::helper::internal::LraViHelper viHelper(bscc, this->_transitionMatrix, aperiodicFactor, this->_markovianStates, this->_exitRates); + return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, this->_exitRates); + } else { + // We assume a DTMC (with deterministic timed states and no instant states) + storm::modelchecker::helper::internal::LraViHelper viHelper(bscc, this->_transitionMatrix, aperiodicFactor); + return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter); + } + } + + template + std::pair, std::vector> SparseDeterministicInfiniteHorizonHelper::buildSspMatrixVector(std::vector const& bsccLraValues, std::vector const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem) { + + // Create SSP Matrix. + // In contrast to the version for nondeterministic models, we eliminate the auxiliary states representing each BSCC on the fly + + // Probability mass that would lead to a BSCC will be considered in the rhs of the equation system + auto sspMatrix = this->_transitionMatrix.getSubmatrix(false, statesNotInComponent, statesNotInComponent, asEquationSystem); + if (asEquationSystem) { + sspMatrix.convertToEquationSystem(); + } + + // Create the SSP right-hand-side + std::vector rhs; + rhs.reserve(sspMatrix.getRowCount()); + for (auto const& state : statesNotInComponent) { + ValueType stateValue = storm::utility::zero(); + for (auto const& transition : this->_transitionMatrix.getRow(state)) { + if (!statesNotInComponent.get(transition.getColumn())) { + // This transition leads to a BSCC! + stateValue += transition.getValue() * bsccLraValues[inputStateToBsccIndexMap[transition.getColumn()]]; + } + } + rhs.push_back(std::move(stateValue)); + } + + return std::make_pair(std::move(sspMatrix), std::move(rhs)); + } + + template + std::vector SparseDeterministicInfiniteHorizonHelper::buildAndSolveSsp(Environment const& env, std::vector const& componentLraValues) { + STORM_LOG_ASSERT(this->_longRunComponentDecomposition != nullptr, "Decomposition not computed, yet."); + + // For fast transition rewriting, we build a mapping from the input state indices to the state indices of a new transition matrix + // which redirects all transitions leading to a former BSCC state to a new (imaginary) auxiliary state. + // Each auxiliary state gets assigned the value of that BSCC and we compute expected rewards (aka stochastic shortest path, SSP) on that new system. + // For efficiency reasons, we actually build the system where the auxiliary states are already eliminated. + + // First gather the states that are part of a component + // and create a mapping from states that lie in a component to the corresponding component index. + storm::storage::BitVector statesInComponents(this->_transitionMatrix.getRowGroupCount()); + std::vector stateIndexMap(this->_transitionMatrix.getRowGroupCount(), std::numeric_limits::max()); + for (uint64_t currentComponentIndex = 0; currentComponentIndex < this->_longRunComponentDecomposition->size(); ++currentComponentIndex) { + for (auto const& element : (*this->_longRunComponentDecomposition)[currentComponentIndex]) { + uint64_t state = internal::getComponentElementState(element); + statesInComponents.set(state); + stateIndexMap[state] = currentComponentIndex; + } + } + // Map the non-component states to their index in the SSP. Note that the order of these states will be preserved. + uint64_t numberOfNonComponentStates = 0; + storm::storage::BitVector statesNotInComponent = ~statesInComponents; + for (auto const& nonComponentState : statesNotInComponent) { + stateIndexMap[nonComponentState] = numberOfNonComponentStates; + ++numberOfNonComponentStates; + } + + // The next step is to create the equation system solving the SSP (unless the whole system consists of BSCCs) + std::vector sspValues; + if (numberOfNonComponentStates > 0) { + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + bool isEqSysFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + auto sspMatrixVector = buildSspMatrixVector(componentLraValues, stateIndexMap, statesNotInComponent, isEqSysFormat); + std::unique_ptr> solver = linearEquationSolverFactory.create(env, sspMatrixVector.first); + auto lowerUpperBounds = std::minmax_element(componentLraValues.begin(), componentLraValues.end()); + solver->setBounds(*lowerUpperBounds.first, *lowerUpperBounds.second); + // Check solver requirements + auto requirements = solver->getRequirements(env); + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + sspValues.assign(sspMatrixVector.first.getRowCount(), (*lowerUpperBounds.first + *lowerUpperBounds.second) / storm::utility::convertNumber(2)); + solver->solveEquations(env, sspValues, sspMatrixVector.second); + } + + // Prepare result vector. + std::vector result(this->_transitionMatrix.getRowGroupCount()); + for (uint64_t state = 0; state < stateIndexMap.size(); ++state) { + if (statesNotInComponent.get(state)) { + result[state] = sspValues[stateIndexMap[state]]; + } else { + result[state] = componentLraValues[stateIndexMap[state]]; + } + } + return result; + } + + template class SparseDeterministicInfiniteHorizonHelper; + template class SparseDeterministicInfiniteHorizonHelper; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h new file mode 100644 index 000000000..5c22cc1af --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h @@ -0,0 +1,65 @@ +#pragma once +#include "storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.h" + + +namespace storm { + + namespace modelchecker { + namespace helper { + + /*! + * Helper class for model checking queries that depend on the long run behavior of the (nondeterministic) system. + * @tparam ValueType the type a value can have + */ + template + class SparseDeterministicInfiniteHorizonHelper : public SparseInfiniteHorizonHelper { + + public: + /*! + * Function mapping from indices to values + */ + typedef typename SparseInfiniteHorizonHelper::ValueGetter ValueGetter; + + /*! + * Initializes the helper for a discrete time model (i.e. DTMC) + */ + SparseDeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix); + + /*! + * Initializes the helper for a continuous time model (i.e. CTMC) + * @note The transition matrix shall be probabilistic (i.e. the rows sum up to one) + */ + SparseDeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates); + + /*! + * @param stateValuesGetter a function returning a value for a given state index + * @param actionValuesGetter a function returning a value for a given (global) choice index + * @return the (unique) optimal LRA value for the given component. + * @post if scheduler production is enabled and Nondeterministic is true, getProducedOptimalChoices() contains choices for the states of the given component which yield the returned LRA value. Choices for states outside of the component are not affected. + */ + virtual ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& component) override; + + protected: + + virtual void createDecomposition() override; + + std::pair computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + + /*! + * As computeLraForMec but uses value iteration as a solution method (independent of what is set in env) + */ + ValueType computeLraForBsccVi(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + + std::pair, std::vector> buildSspMatrixVector(std::vector const& bsccLraValues, std::vector const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem); + + /*! + * @return Lra values for each state + */ + virtual std::vector buildAndSolveSsp(Environment const& env, std::vector const& mecLraValues) override; + + }; + + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.cpp new file mode 100644 index 000000000..8eab0bca0 --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.cpp @@ -0,0 +1,161 @@ +#include "SparseInfiniteHorizonHelper.h" + +#include "storm/modelchecker/helper/infinitehorizon/internal/ComponentUtility.h" +#include "storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h" + +#include "storm/models/sparse/StandardRewardModel.h" + +#include "storm/storage/SparseMatrix.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" +#include "storm/storage/StronglyConnectedComponentDecomposition.h" + +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/Multiplier.h" +#include "storm/solver/LpSolver.h" + +#include "storm/utility/SignalHandler.h" +#include "storm/utility/solver.h" +#include "storm/utility/vector.h" + +#include "storm/environment/solver/LongRunAverageSolverEnvironment.h" +#include "storm/environment/solver/MinMaxSolverEnvironment.h" + +#include "storm/exceptions/UnmetRequirementException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + SparseInfiniteHorizonHelper::SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix) : _transitionMatrix(transitionMatrix), _markovianStates(nullptr), _exitRates(nullptr), _backwardTransitions(nullptr), _longRunComponentDecomposition(nullptr) { + // Intentionally left empty. + } + + template + SparseInfiniteHorizonHelper::SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates) : _transitionMatrix(transitionMatrix), _markovianStates(&markovianStates), _exitRates(&exitRates), _backwardTransitions(nullptr), _longRunComponentDecomposition(nullptr) { + // Intentionally left empty. + } + + template + SparseInfiniteHorizonHelper::SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates) : _transitionMatrix(transitionMatrix), _markovianStates(nullptr), _exitRates(&exitRates), _backwardTransitions(nullptr), _longRunComponentDecomposition(nullptr) { + // Intentionally left empty. + } + + template + void SparseInfiniteHorizonHelper::provideBackwardTransitions(storm::storage::SparseMatrix const& backwardTransitions) { + STORM_LOG_WARN_COND(_backwardTransitions == nullptr, "Backwards transitions were provided but they were already computed or provided before."); + _backwardTransitions = &backwardTransitions; + } + + template + void SparseInfiniteHorizonHelper::provideLongRunComponentDecomposition(storm::storage::Decomposition const& decomposition) { + STORM_LOG_WARN_COND(_longRunComponentDecomposition == nullptr, "Long Run Component Decomposition was provided but it was already computed or provided before."); + _longRunComponentDecomposition = &decomposition; + } + + template + std::vector SparseInfiniteHorizonHelper::computeLongRunAverageProbabilities(Environment const& env, storm::storage::BitVector const& psiStates) { + return computeLongRunAverageValues(env, + [&psiStates] (uint64_t stateIndex) { return psiStates.get(stateIndex) ? storm::utility::one() : storm::utility::zero(); }, + [] (uint64_t) { return storm::utility::zero(); } + ); + } + + template + std::vector SparseInfiniteHorizonHelper::computeLongRunAverageRewards(Environment const& env, storm::models::sparse::StandardRewardModel const& rewardModel) { + ValueGetter stateRewardsGetter; + if (rewardModel.hasStateRewards()) { + stateRewardsGetter = [&rewardModel] (uint64_t stateIndex) { return rewardModel.getStateReward(stateIndex); }; + } else { + stateRewardsGetter = [] (uint64_t) { return storm::utility::zero(); }; + } + ValueGetter actionRewardsGetter; + if (rewardModel.hasStateActionRewards() || rewardModel.hasTransitionRewards()) { + if (rewardModel.hasTransitionRewards()) { + actionRewardsGetter = [&] (uint64_t globalChoiceIndex) { return rewardModel.getStateActionAndTransitionReward(globalChoiceIndex, this->_transitionMatrix); }; + } else { + actionRewardsGetter = [&] (uint64_t globalChoiceIndex) { return rewardModel.getStateActionReward(globalChoiceIndex); }; + } + } else { + stateRewardsGetter = [] (uint64_t) { return storm::utility::zero(); }; + } + + return computeLongRunAverageValues(env, stateRewardsGetter, actionRewardsGetter); + } + + template + std::vector SparseInfiniteHorizonHelper::computeLongRunAverageValues(Environment const& env, std::vector const* stateValues, std::vector const* actionValues) { + ValueGetter stateValuesGetter; + if (stateValues) { + stateValuesGetter = [&stateValues] (uint64_t stateIndex) { return (*stateValues)[stateIndex]; }; + } else { + stateValuesGetter = [] (uint64_t) { return storm::utility::zero(); }; + } + ValueGetter actionValuesGetter; + if (actionValues) { + actionValuesGetter = [&actionValues] (uint64_t globalChoiceIndex) { return (*actionValues)[globalChoiceIndex]; }; + } else { + actionValuesGetter = [] (uint64_t) { return storm::utility::zero(); }; + } + + return computeLongRunAverageValues(env, stateValuesGetter, actionValuesGetter); + + } + + + template + std::vector SparseInfiniteHorizonHelper::computeLongRunAverageValues(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter) { + // We will compute the long run average value for each MEC individually and then set-up an Equation system to compute the value also at non-mec states. + // For a description of this approach see, e.g., Guck et al.: Modelling and Analysis of Markov Reward Automata (ATVA'14), https://doi.org/10.1007/978-3-319-11936-6_13 + + // Prepare an environment for the underlying solvers + auto underlyingSolverEnvironment = env; + if (env.solver().isForceSoundness()) { + // For sound computations, the error in the MECS plus the error in the remaining system should not exceed the user defined precsion. + storm::RationalNumber newPrecision = env.solver().lra().getPrecision() / storm::utility::convertNumber(2); + underlyingSolverEnvironment.solver().minMax().setPrecision(newPrecision); + underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion()); + underlyingSolverEnvironment.solver().setLinearEquationSolverPrecision(newPrecision, env.solver().lra().getRelativeTerminationCriterion()); + underlyingSolverEnvironment.solver().lra().setPrecision(newPrecision); + } + + // If requested, allocate memory for the choices made + if (Nondeterministic && this->isProduceSchedulerSet()) { + if (!_producedOptimalChoices.is_initialized()) { + _producedOptimalChoices.emplace(); + } + _producedOptimalChoices->resize(_transitionMatrix.getRowGroupCount()); + } + STORM_LOG_ASSERT(Nondeterministic || !this->isProduceSchedulerSet(), "Scheduler production enabled for deterministic model."); + + // Decompose the model to their bottom components (MECS or BSCCS) + createDecomposition(); + + // Compute the long-run average for all components in isolation. + std::vector componentLraValues; + componentLraValues.reserve(_longRunComponentDecomposition->size()); + for (auto const& c : *_longRunComponentDecomposition) { + componentLraValues.push_back(computeLraForComponent(underlyingSolverEnvironment, stateRewardsGetter, actionRewardsGetter, c)); + } + + // Solve the resulting SSP where end components are collapsed into single auxiliary states + return buildAndSolveSsp(underlyingSolverEnvironment, componentLraValues); + } + + template + bool SparseInfiniteHorizonHelper::isContinuousTime() const { + STORM_LOG_ASSERT((_markovianStates == nullptr) || (_exitRates != nullptr), "Inconsistent information given: Have Markovian states but no exit rates." ); + return _exitRates != nullptr; + } + + template class SparseInfiniteHorizonHelper; + template class SparseInfiniteHorizonHelper; + template class SparseInfiniteHorizonHelper; + + template class SparseInfiniteHorizonHelper; + template class SparseInfiniteHorizonHelper; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.h b/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.h new file mode 100644 index 000000000..b8ebcef0b --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.h @@ -0,0 +1,140 @@ +#pragma once +#include "storm/modelchecker/helper/SingleValueModelCheckerHelper.h" + +#include "storm/storage/MaximalEndComponent.h" +#include "storm/storage/StronglyConnectedComponent.h" +#include "storm/storage/Decomposition.h" +#include "storm/storage/SparseMatrix.h" + +namespace storm { + class Environment; + + namespace models { + namespace sparse { + template class StandardRewardModel; + } + } + + namespace modelchecker { + namespace helper { + + /*! + * Helper class for model checking queries that depend on the long run behavior of the (nondeterministic) system. + * @tparam ValueType the type a value can have + * @tparam Nondeterministic true if there is nondeterminism in the Model (MDP or MA) + */ + template + class SparseInfiniteHorizonHelper : public SingleValueModelCheckerHelper { + + public: + + /*! + * The type of a component in which the system resides in the long run (BSCC for deterministic models, MEC for nondeterministic models) + */ + using LongRunComponentType = typename std::conditional::type; + + /*! + * Function mapping from indices to values + */ + typedef std::function ValueGetter; + + /*! + * Initializes the helper for a discrete time (i.e. DTMC, MDP) + */ + SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix); + + /*! + * Initializes the helper for continuous time (i.e. MA) + */ + SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates); + + /*! + * Initializes the helper for continuous time (i.e. CTMC) + */ + SparseInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates); + + /*! + * Provides backward transitions that can be used during the computation. + * Providing them is optional. If they are not provided, they will be computed internally + * Be aware that this class does not take ownership, i.e. the caller has to make sure that the reference to the backwardstransitions remains valid. + */ + void provideBackwardTransitions(storm::storage::SparseMatrix const& backwardsTransitions); + + /*! + * Provides the decomposition into long run components (BSCCs/MECs) that can be used during the computation. + * Providing the decomposition is optional. If it is not provided, they will be computed internally. + * Be aware that this class does not take ownership, i.e. the caller has to make sure that the reference to the decomposition remains valid. + */ + void provideLongRunComponentDecomposition(storm::storage::Decomposition const& decomposition); + + /*! + * Computes the long run average probabilities, i.e., the fraction of the time we are in a psiState + * @return a value for each state + */ + std::vector computeLongRunAverageProbabilities(Environment const& env, storm::storage::BitVector const& psiStates); + + /*! + * Computes the long run average rewards, i.e., the average reward collected per time unit + * @return a value for each state + */ + std::vector computeLongRunAverageRewards(Environment const& env, storm::models::sparse::StandardRewardModel const& rewardModel); + + /*! + * Computes the long run average value given the provided state and action-based rewards. + * @param stateValues a vector containing a value for every state + * @param actionValues a vector containing a value for every choice + * @return a value for each state + */ + std::vector computeLongRunAverageValues(Environment const& env, std::vector const* stateValues = nullptr, std::vector const* actionValues = nullptr); + + /*! + * Computes the long run average value given the provided state and action based rewards + * @param stateValuesGetter a function returning a value for a given state index + * @param actionValuesGetter a function returning a value for a given (global) choice index + * @return a value for each state + */ + std::vector computeLongRunAverageValues(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter); + + /*! + * @param stateValuesGetter a function returning a value for a given state index + * @param actionValuesGetter a function returning a value for a given (global) choice index + * @return the (unique) optimal LRA value for the given component. + * @post if scheduler production is enabled and Nondeterministic is true, getProducedOptimalChoices() contains choices for the states of the given component which yield the returned LRA value. Choices for states outside of the component are not affected. + */ + virtual ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& component) = 0; + + protected: + + /*! + * @return true iff this is a computation on a continuous time model (i.e. CTMC, MA) + */ + bool isContinuousTime() const; + + /*! + * @post _longRunComponentDecomposition points to a decomposition of the long run components (MECs, BSCCs) + */ + virtual void createDecomposition() = 0; + + /*! + * @pre if scheduler production is enabled and Nondeterministic is true, a choice for each state within a component must be set such that the choices yield optimal values w.r.t. the individual components. + * @return Lra values for each state + * @post if scheduler production is enabled and Nondeterministic is true, getProducedOptimalChoices() contains choices for all input model states which yield the returned LRA values. + */ + virtual std::vector buildAndSolveSsp(Environment const& env, std::vector const& mecLraValues) = 0; + + storm::storage::SparseMatrix const& _transitionMatrix; + storm::storage::BitVector const* _markovianStates; + std::vector const* _exitRates; + + storm::storage::SparseMatrix const* _backwardTransitions; + storm::storage::Decomposition const* _longRunComponentDecomposition; + std::unique_ptr> _computedBackwardTransitions; + std::unique_ptr> _computedLongRunComponentDecomposition; + + boost::optional> _producedOptimalChoices; + }; + + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp index 5cb86d83e..4dcd8088f 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp @@ -3,18 +3,14 @@ #include "storm/modelchecker/helper/infinitehorizon/internal/ComponentUtility.h" #include "storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h" -#include "storm/models/sparse/StandardRewardModel.h" - #include "storm/storage/SparseMatrix.h" #include "storm/storage/MaximalEndComponentDecomposition.h" -#include "storm/storage/StronglyConnectedComponentDecomposition.h" +#include "storm/storage/Scheduler.h" #include "storm/solver/MinMaxLinearEquationSolver.h" -#include "storm/solver/LinearEquationSolver.h" #include "storm/solver/Multiplier.h" #include "storm/solver/LpSolver.h" -#include "storm/utility/SignalHandler.h" #include "storm/utility/solver.h" #include "storm/utility/vector.h" @@ -27,146 +23,32 @@ namespace storm { namespace modelchecker { namespace helper { - template - SparseNondeterministicInfiniteHorizonHelper::SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix) : _transitionMatrix(transitionMatrix), _backwardTransitions(nullptr), _longRunComponentDecomposition(nullptr), _markovianStates(nullptr), _exitRates(nullptr) { + template + SparseNondeterministicInfiniteHorizonHelper::SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix) : SparseInfiniteHorizonHelper(transitionMatrix) { // Intentionally left empty. } - template - SparseNondeterministicInfiniteHorizonHelper::SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates) : _transitionMatrix(transitionMatrix), _backwardTransitions(nullptr), _longRunComponentDecomposition(nullptr), _markovianStates(&markovianStates), _exitRates(&exitRates) { + template + SparseNondeterministicInfiniteHorizonHelper::SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates) : SparseInfiniteHorizonHelper(transitionMatrix, markovianStates, exitRates) { // Intentionally left empty. } - template - void SparseNondeterministicInfiniteHorizonHelper::provideBackwardTransitions(storm::storage::SparseMatrix const& backwardTransitions) { - STORM_LOG_WARN_COND(_backwardTransitions == nullptr, "Backwards transitions were provided but they were already computed or provided before."); - _backwardTransitions = &backwardTransitions; - } - - template - void SparseNondeterministicInfiniteHorizonHelper::provideLongRunComponentDecomposition(storm::storage::Decomposition const& decomposition) { - STORM_LOG_WARN_COND(_longRunComponentDecomposition == nullptr, "Long Run Component Decomposition was provided but it was already computed or provided before."); - _longRunComponentDecomposition = &decomposition; - } - - template - std::vector SparseNondeterministicInfiniteHorizonHelper::computeLongRunAverageProbabilities(Environment const& env, storm::storage::BitVector const& psiStates) { - return computeLongRunAverageValues(env, - [&psiStates] (uint64_t stateIndex) { return psiStates.get(stateIndex) ? storm::utility::one() : storm::utility::zero(); }, - [] (uint64_t) { return storm::utility::zero(); } - ); - } - - template - std::vector SparseNondeterministicInfiniteHorizonHelper::computeLongRunAverageRewards(Environment const& env, storm::models::sparse::StandardRewardModel const& rewardModel) { - std::function stateRewardsGetter; - if (rewardModel.hasStateRewards()) { - stateRewardsGetter = [&rewardModel] (uint64_t stateIndex) { return rewardModel.getStateReward(stateIndex); }; - } else { - stateRewardsGetter = [] (uint64_t) { return storm::utility::zero(); }; - } - std::function actionRewardsGetter; - if (rewardModel.hasStateActionRewards() || rewardModel.hasTransitionRewards()) { - if (rewardModel.hasTransitionRewards()) { - actionRewardsGetter = [&] (uint64_t globalChoiceIndex) { return rewardModel.getStateActionAndTransitionReward(globalChoiceIndex, this->_transitionMatrix); }; - } else { - actionRewardsGetter = [&] (uint64_t globalChoiceIndex) { return rewardModel.getStateActionReward(globalChoiceIndex); }; - } - } else { - stateRewardsGetter = [] (uint64_t) { return storm::utility::zero(); }; - } - - return computeLongRunAverageValues(env, stateRewardsGetter, actionRewardsGetter); - } - - template - std::vector SparseNondeterministicInfiniteHorizonHelper::computeLongRunAverageValues(Environment const& env, std::vector const* stateValues, std::vector const* actionValues) { - std::function stateValuesGetter; - if (stateValues) { - stateValuesGetter = [&stateValues] (uint64_t stateIndex) { return (*stateValues)[stateIndex]; }; - } else { - stateValuesGetter = [] (uint64_t) { return storm::utility::zero(); }; - } - std::function actionValuesGetter; - if (actionValues) { - actionValuesGetter = [&actionValues] (uint64_t globalChoiceIndex) { return (*actionValues)[globalChoiceIndex]; }; - } else { - actionValuesGetter = [] (uint64_t) { return storm::utility::zero(); }; - } - - return computeLongRunAverageValues(env, stateValuesGetter, actionValuesGetter); - - } - - template - std::vector SparseNondeterministicInfiniteHorizonHelper::computeLongRunAverageValues(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter) { - // We will compute the long run average value for each MEC individually and then set-up a MinMax Equation system to compute the value also at non-mec states. - // For a description of this approach see, e.g., Guck et al.: Modelling and Analysis of Markov Reward Automata (ATVA'14), https://doi.org/10.1007/978-3-319-11936-6_13 - - // Prepare an environment for the underlying solvers - auto underlyingSolverEnvironment = env; - if (env.solver().isForceSoundness()) { - // For sound computations, the error in the MECS plus the error in the remaining system should not exceed the user defined precsion. - storm::RationalNumber newPrecision = env.solver().lra().getPrecision() / storm::utility::convertNumber(2); - underlyingSolverEnvironment.solver().minMax().setPrecision(newPrecision); - underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion()); - underlyingSolverEnvironment.solver().setLinearEquationSolverPrecision(newPrecision, env.solver().lra().getRelativeTerminationCriterion()); - underlyingSolverEnvironment.solver().lra().setPrecision(newPrecision); - } - - // If requested, allocate memory for the choices made - if (Nondeterministic && this->isProduceSchedulerSet()) { - if (!_producedOptimalChoices.is_initialized()) { - _producedOptimalChoices.emplace(); - } - _producedOptimalChoices->resize(_transitionMatrix.getRowGroupCount()); - } - STORM_LOG_ASSERT(Nondeterministic || !this->isProduceSchedulerSet(), "Scheduler production enabled for deterministic model."); - - // Start by decomposing the Model into its MECs. - if (_longRunComponentDecomposition == nullptr) { - // The decomposition has not been provided or computed, yet. - if (Nondeterministic) { - if (_backwardTransitions == nullptr) { - _computedBackwardTransitions = std::make_unique(_transitionMatrix.transpose(true)); - _backwardTransitions = _computedBackwardTransitions.get(); - } - _computedLongRunComponentDecomposition = std::make_unique(_transitionMatrix, *_backwardTransitions); - } else { - _computedLongRunComponentDecomposition = std::make_unique(_transitionMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs()); - } - _longRunComponentDecomposition = _computedLongRunComponentDecomposition.get(); - } - - // Compute the long-run average for all components in isolation. - std::vector componentLraValues; - mecLraValues.reserve(_longRunComponentDecomposition->size()); - for (auto const& c : *_longRunComponentDecomposition) { - componentLraValues.push_back(computeLraForComponent(underlyingSolverEnvironment, stateRewardsGetter, actionRewardsGetter, c)); - } - - // Solve the resulting SSP where end components are collapsed into single auxiliary states - return buildAndSolveSsp(underlyingSolverEnvironment, componentLraValues); - } - - template - std::vector const& SparseNondeterministicInfiniteHorizonHelper::getProducedOptimalChoices() const { - STORM_LOG_WARN_COND(Nondeterministic, "Getting optimal choices for deterministic model."); + template + std::vector const& SparseNondeterministicInfiniteHorizonHelper::getProducedOptimalChoices() const { STORM_LOG_ASSERT(this->isProduceSchedulerSet(), "Trying to get the produced optimal choices although no scheduler was requested."); - STORM_LOG_ASSERT(_producedOptimalChoices.is_initialized(), "Trying to get the produced optimal choices but none were available. Was there a computation call before?"); - return _producedOptimalChoices.get(); + 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& SparseNondeterministicInfiniteHorizonHelper::getProducedOptimalChoices() { - STORM_LOG_WARN_COND(Nondeterministic, "Getting optimal choices for deterministic model."); + template + std::vector& SparseNondeterministicInfiniteHorizonHelper::getProducedOptimalChoices() { STORM_LOG_ASSERT(this->isProduceSchedulerSet(), "Trying to get the produced optimal choices although no scheduler was requested."); - STORM_LOG_ASSERT(_producedOptimalChoices.is_initialized(), "Trying to get the produced optimal choices but none were available. Was there a computation call before?"); - return _producedOptimalChoices.get(); + 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 SparseNondeterministicInfiniteHorizonHelper::extractScheduler() const { + template + storm::storage::Scheduler SparseNondeterministicInfiniteHorizonHelper::extractScheduler() const { auto const& optimalChoices = getProducedOptimalChoices(); storm::storage::Scheduler scheduler(optimalChoices.size()); for (uint64_t state = 0; state < optimalChoices.size(); ++state) { @@ -175,44 +57,32 @@ namespace storm { return scheduler; } - template - bool SparseNondeterministicInfiniteHorizonHelper::isContinuousTime() const { - STORM_LOG_ASSERT((_markovianStates == nullptr) || (_exitRates != nullptr), "Inconsistent information given: Have Markovian states but no exit rates." ); - return _exitRates != nullptr; - } - - template - template < typename = typename std::enable_if< !Nondeterministic >::type > - ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, LongRunComponentType const& component) { - // For deterministic models, we compute the LRA for a BSCC - - STORM_LOG_ASSERT(!this->isProduceSchedulerSet(), "Scheduler production enabled for deterministic model."); - - auto trivialResult = computeLraForTrivialComponent(env, stateReardsGetter, actionRewardsGetter, component); - if (trivialResult.first) { - return trivialResult.second; + template + void SparseNondeterministicInfiniteHorizonHelper::createDecomposition() { + if (this->_longRunComponentDecomposition == nullptr) { + // The decomposition has not been provided or computed, yet. + if (this->_backwardTransitions == nullptr) { + this->_computedBackwardTransitions = std::make_unique>(this->_transitionMatrix.transpose(true)); + this->_backwardTransitions = this->_computedBackwardTransitions.get(); + } + this->_computedLongRunComponentDecomposition = std::make_unique>(this->_transitionMatrix, *this->_backwardTransitions); + this->_longRunComponentDecomposition = this->_computedLongRunComponentDecomposition.get(); } - - // Solve nontrivial BSCC with the method specified in the settings - - // TODO - } - - template - template < typename = typename std::enable_if< Nondeterministic >::type > - ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, LongRunComponentType const& component) { + + template + ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::MaximalEndComponent const& component) { // For models with potential nondeterminisim, we compute the LRA for a maximal end component (MEC) // Allocate memory for the nondeterministic choices. if (this->isProduceSchedulerSet()) { - if (!_producedOptimalChoices.is_initialized()) { - _producedOptimalChoices.emplace(); + if (!this->_producedOptimalChoices.is_initialized()) { + this->_producedOptimalChoices.emplace(); } - _producedOptimalChoices->resize(_transitionMatrix.getRowGroupCount()); + this->_producedOptimalChoices->resize(this->_transitionMatrix.getRowGroupCount()); } - auto trivialResult = computeLraForTrivialComponent(env, stateReardsGetter, actionRewardsGetter, component); + auto trivialResult = this->computeLraForTrivialMec(env, stateRewardsGetter, actionRewardsGetter, component); if (trivialResult.first) { return trivialResult.second; } @@ -228,28 +98,28 @@ namespace storm { } STORM_LOG_ERROR_COND(!this->isProduceSchedulerSet() || method == storm::solver::LraMethod::ValueIteration, "Scheduler generation not supported for the chosen LRA method. Try value-iteration."); if (method == storm::solver::LraMethod::LinearProgramming) { - return computeLraForMecLp(env, stateRewardsGetter, actionRewardsGetter, mec); + return computeLraForMecLp(env, stateRewardsGetter, actionRewardsGetter, component); } else if (method == storm::solver::LraMethod::ValueIteration) { - return computeLraForMecVi(env, stateRewardsGetter, actionRewardsGetter, mec); + return computeLraForMecVi(env, stateRewardsGetter, actionRewardsGetter, component); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } } - - template - std::pair SparseNondeterministicInfiniteHorizonHelper::computeLraForTrivialComponent(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, LongRunComponentType const& component) { + + template + std::pair SparseNondeterministicInfiniteHorizonHelper::computeLraForTrivialMec(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::MaximalEndComponent const& component) { // If the component only consists of a single state, we compute the LRA value directly if (component.size() == 1) { - auto element const& = *component.begin(); + auto const& element = *component.begin(); uint64_t state = internal::getComponentElementState(element); - auto choiceIt = internal::getComponentChoicesBegin(element); - if (Nondeterministic && !isContinuousTime()) { + auto choiceIt = internal::getComponentElementChoicesBegin(element); + if (!this->isContinuousTime()) { // This is an MDP. // Find the choice with the highest/lowest reward ValueType bestValue = actionRewardsGetter(*choiceIt); uint64_t bestChoice = *choiceIt; - for (++choiceIt; choiceIt != getComponentChoicesEnd(element); ++choiceIt) { + for (++choiceIt; choiceIt != internal::getComponentElementChoicesEnd(element); ++choiceIt) { ValueType currentValue = actionRewardsGetter(*choiceIt); if ((this->minimize() && currentValue < bestValue) || (this->maximize() && currentValue > bestValue)) { bestValue = std::move(currentValue); @@ -257,65 +127,50 @@ namespace storm { } } if (this->isProduceSchedulerSet()) { - _producedOptimalChoices.get()[state] = bestChoice - _transitionMatrix.getRowGroupIndices()[state]; + this->_producedOptimalChoices.get()[state] = bestChoice - this->_transitionMatrix.getRowGroupIndices()[state]; } bestValue += stateRewardsGetter(state); return {true, bestValue}; } else { // In a Markov Automaton, singleton components have to consist of a Markovian state because of the non-Zenoness assumption. Then, there is just one possible choice. - STORM_LOG_THROW(!Nondeterministic || (_markovianStates != nullptr && _markovianStates->get(state)), storm::exceptions::InvalidOperationException, "Markov Automaton has Zeno behavior. Computation of Long Run Average values not supported."); + STORM_LOG_ASSERT(this->_markovianStates != nullptr, "Nondeterministic continuous time model without Markovian states... Is this a not a Markov Automaton?"); + STORM_LOG_THROW(this->_markovianStates->get(state), storm::exceptions::InvalidOperationException, "Markov Automaton has Zeno behavior. Computation of Long Run Average values not supported."); STORM_LOG_ASSERT(internal::getComponentElementChoiceCount(element) == 1, "Markovian state has Nondeterministic behavior."); - if (Nondeterministic && this->isProduceSchedulerSet()) { - _producedOptimalChoices.get()[state] = 0; + if (this->isProduceSchedulerSet()) { + this->_producedOptimalChoices.get()[state] = 0; } - ValueType result = stateRewardsGetter(state) + (isContinuousTime() ? (*_exitRates)[state] * actionRewardsGetter(*choiceIt) : actionRewardsGetter(*choiceIt)); + ValueType result = stateRewardsGetter(state) + (this->isContinuousTime() ? (*this->_exitRates)[state] * actionRewardsGetter(*choiceIt) : actionRewardsGetter(*choiceIt)); return {true, result}; } - } else if (!Nondeterministic) { - // For deterministic models, we can also easily catch the case where all values are the same - bool first = true; - ValueType val = storm::utility::zero(); - for (auto const& element : component) { - auto state = getComponentElementState(element); - STORM_LOG_ASSERT(state == *getComponentChoicesBegin(element), "Unexpected choice index at state " << state << " of deterministic model."); - ValueType curr = stateRewardsGetter(state) + (isContinuousTime() ? (*_exitRates)[state] * actionRewardsGetter(state) : actionRewardsGetter(state)); - if (first) { - first = false; - } else if (val != curr) { - return {false, storm::utility::zero()}; - } - } - // All values are the same - return {true, val}; - } else { - return {false, storm::utility::zero()}; } + return {false, storm::utility::zero()}; } - template - ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForMecVi(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) { + + template + ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForMecVi(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) { // Collect some parameters of the computation ValueType aperiodicFactor = storm::utility::convertNumber(env.solver().lra().getAperiodicFactor()); std::vector* optimalChoices = nullptr; if (this->isProduceSchedulerSet()) { - optimalChoices = &_producedOptimalChoices.get(); + optimalChoices = &this->_producedOptimalChoices.get(); } // Now create a helper and perform the algorithm - if (isContinuousTime()) { + if (this->isContinuousTime()) { // We assume a Markov Automaton (with deterministic timed states and nondeterministic instant states) - storm::modelchecker::helper::internal::LraViHelper viHelper(mec, _transitionMatrix, aperiodicFactor, _markovianStates, _exitRates); - return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, _exitRates, &this->getOptimizationDirection(), optimalChoices); + storm::modelchecker::helper::internal::LraViHelper viHelper(mec, this->_transitionMatrix, aperiodicFactor, this->_markovianStates, this->_exitRates); + return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, this->_exitRates, &this->getOptimizationDirection(), optimalChoices); } else { // We assume an MDP (with nondeterministic timed states and no instant states) - storm::modelchecker::helper::internal::LraViHelper viHelper(mec, _transitionMatrix, aperiodicFactor); + storm::modelchecker::helper::internal::LraViHelper viHelper(mec, this->_transitionMatrix, aperiodicFactor); return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, nullptr, &this->getOptimizationDirection(), optimalChoices); } } - template - ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForMecLp(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) { + template + ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForMecLp(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) { // Create an LP solver auto solver = storm::utility::solver::LpSolverFactory().create("LRA for MEC"); @@ -336,18 +191,18 @@ namespace storm { // Add constraints. for (auto const& stateChoicesPair : mec) { uint_fast64_t state = stateChoicesPair.first; - bool stateIsMarkovian = _markovianStates && _markovianStates->get(state); + bool stateIsMarkovian = this->_markovianStates && this->_markovianStates->get(state); // Now create a suitable constraint for each choice // x_s {≤, ≥} -k/rate(s) + sum_s' P(s,act,s') * x_s' + (value(s)/rate(s) + value(s,act)) for (auto choice : stateChoicesPair.second) { std::vector summands; - auto matrixRow = _transitionMatrix.getRow(choice); + auto matrixRow = this->_transitionMatrix.getRow(choice); summands.reserve(matrixRow.getNumberOfEntries() + 2); // add -k/rate(s) (only if s is either a Markovian state or we have an MDP) if (stateIsMarkovian) { - summands.push_back(-(k / solver->getManager().rational((*_exitRates)[state]))); - } else if (!isContinuousTime()) { + summands.push_back(-(k / solver->getManager().rational((*this->_exitRates)[state]))); + } else if (!this->isContinuousTime()) { summands.push_back(-k); } // add sum_s' P(s,act,s') * x_s' @@ -358,8 +213,8 @@ namespace storm { ValueType value; if (stateIsMarkovian) { // divide state reward with exit rate - value = stateRewardsGetter(state) / (*_exitRates)[state] + actionRewardsGetter(choice); - } else if (!isContinuousTime()) { + value = stateRewardsGetter(state) / (*this->_exitRates)[state] + actionRewardsGetter(choice); + } else if (!this->isContinuousTime()) { // in discrete time models no scaling is needed value = stateRewardsGetter(state) + actionRewardsGetter(choice); } else { @@ -388,7 +243,7 @@ namespace storm { * Transitions that don't lead to a Component state are copied (taking a state index mapping into account). */ template - void addSspMatrixChoice(uint64_t const& inputMatrixChoice, storm::storage::SparseMatrix const& inputTransitionMatrix, std::vector const& inputToSspStateMap, uint64_t const& numberOfStatesNotInComponents, uint64_t const& currentSspChoice, storm::storage::SparseMatrixBuilder& sspMatrixBuilder) { + void addSspMatrixChoice(uint64_t const& inputMatrixChoice, storm::storage::SparseMatrix const& inputTransitionMatrix, std::vector const& inputToSspStateMap, uint64_t const& numberOfNonComponentStates, uint64_t const& currentSspChoice, storm::storage::SparseMatrixBuilder& sspMatrixBuilder) { // As there could be multiple transitions to the same MEC, we accumulate them in this map before adding them to the matrix builder. std::map auxiliaryStateToProbabilityMap; @@ -398,7 +253,7 @@ namespace storm { auto const& sspTransitionTarget = inputToSspStateMap[transition.getColumn()]; // Since the auxiliary Component states are appended at the end of the matrix, we can use this check to // decide whether the transition leads to a component state or not - if (sspTransitionTarget < numberOfStatesNotInMecs) { + if (sspTransitionTarget < numberOfNonComponentStates) { // If the target state is not contained in a component, we can copy over the entry. sspMatrixBuilder.addNextValue(currentSspChoice, sspTransitionTarget, transition.getValue()); } else { @@ -420,195 +275,204 @@ namespace storm { } } - template - std::vector SparseNondeterministicInfiniteHorizonHelper::buildAndSolveSsp(Environment const& env, std::vector const& mecLraValues) { - STORM_LOG_ASSERT(_longRunComponentDecomposition != nullptr, "Decomposition not computed, yet."); - - uint64_t numberOfStates = _transitionMatrix.getRowGroupCount(); - - // For fast transition rewriting, we build a mapping from the input state indices to the state indices of a new transition matrix - // which redirects all transitions leading to a former component state to a new auxiliary state. - // There will be one auxiliary state for each component. These states will be appended to the end of the matrix. - - // First gather the states that are part of a component - // and create a mapping from states that lie in a component to the corresponding component index. - storm::storage::BitVector statesInMecs(numberOfStates); - std::vector inputToSspStateMap(numberOfStates, std::numeric_limits::max()); - for (uint64_t currentMecIndex = 0; currentMecIndex < _longRunComponentDecomposition->size(); ++currentMecIndex) { - for (auto const& stateChoicesPair : (*_longRunComponentDecomposition)[currentMecIndex]) { - statesInMecs.set(stateChoicesPair.first); - inputToSspStateMap[stateChoicesPair.first] = currentMecIndex; - } - } - // Now take care of the non-mec states. Note that the order of these states will be preserved. - uint64_t numberOfStatesNotInMecs = 0; - storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; - for (auto const& nonMecState : statesNotContainedInAnyMec) { - inputToSspStateMap[nonMecState] = numberOfStatesNotInMecs; - ++numberOfStatesNotInMecs; - } - // Finalize the mapping for the mec states which now still assigns mec states to to their Mec index. - // To make sure that they point to the auxiliary states (located at the end of the SspMatrix), we need to shift them by the - // number of states that are not in a mec. - for (auto const& mecState : statesInMecs) { - inputToSspStateMap[mecState] += numberOfStatesNotInMecs; - } + template + std::pair, std::vector> SparseNondeterministicInfiniteHorizonHelper::buildSspMatrixVector(std::vector const& mecLraValues, std::vector const& inputToSspStateMap, storm::storage::BitVector const& statesNotInComponent, uint64_t numberOfNonComponentStates, std::vector>* sspComponentExitChoicesToOriginalMap) { - // For scheduler extraction, we will need to create a mapping between choices at the auxiliary states and the - // corresponding choices in the original model. - std::vector> sspMecExitChoicesToOriginalMap; + auto const& choiceIndices = this->_transitionMatrix.getRowGroupIndices(); - // The next step is to create the SSP matrix and the right-hand side of the SSP. std::vector rhs; - uint64_t numberOfSspStates = numberOfStatesNotInMecs + _longRunComponentDecomposition->size(); - typename storm::storage::SparseMatrixBuilder sspMatrixBuilder(0, numberOfSspStates , 0, false, true, numberOfSspStates); - // If the source state of a transition is not contained in any MEC, we copy its choices (and perform the necessary modifications). + uint64_t numberOfSspStates = numberOfNonComponentStates + this->_longRunComponentDecomposition->size(); + storm::storage::SparseMatrixBuilder sspMatrixBuilder(0, numberOfSspStates , 0, true, true, numberOfSspStates); + // If the source state of a transition is not contained in any component, we copy its choices (and perform the necessary modifications). uint64_t currentSspChoice = 0; - for (auto const& nonMecState : statesNotContainedInAnyMec) { + for (auto const& nonComponentState : statesNotInComponent) { sspMatrixBuilder.newRowGroup(currentSspChoice); - - for (uint64_t choice = nondeterministicChoiceIndices[nonMecState]; choice < nondeterministicChoiceIndices[nonMecState + 1]; ++choice, ++currentSspChoice) { + for (uint64_t choice = choiceIndices[nonComponentState]; choice < choiceIndices[nonComponentState + 1]; ++choice, ++currentSspChoice) { rhs.push_back(storm::utility::zero()); - addSspMatrixChoice(choice, _transitionMatrix, inputToSspStateMap, numberOfStatesNotInMecs, currentSspChoice, sspMatrixBuilder); + addSspMatrixChoice(choice, this->_transitionMatrix, inputToSspStateMap, numberOfNonComponentStates, currentSspChoice, sspMatrixBuilder); } } - // Now we construct the choices for the auxiliary states which reflect former MEC states. - for (uint64_t mecIndex = 0; mecIndex < _longRunComponentDecomposition->size(); ++mecIndex) { - storm::storage::MaximalEndComponent const& mec = (*_longRunComponentDecomposition)[mecIndex]; + // Now we construct the choices for the auxiliary states which reflect former Component states. + for (uint64_t componentIndex = 0; componentIndex < this->_longRunComponentDecomposition->size(); ++componentIndex) { + auto const& component = (*this->_longRunComponentDecomposition)[componentIndex]; sspMatrixBuilder.newRowGroup(currentSspChoice); - for (auto const& stateChoicesPair : mec) { - uint64_t const& mecState = stateChoicesPair.first; - auto const& choicesInMec = stateChoicesPair.second; - for (uint64_t choice = nondeterministicChoiceIndices[mecState]; choice < nondeterministicChoiceIndices[mecState + 1]; ++choice) { - // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state. - if (choicesInMec.find(choice) == choicesInMec.end()) { + // For nondeterministic models it might still be that we leave the component again. This needs to be reflected in the SSP + // by adding the "exiting" choices of the MEC to the axiliary states + for (auto const& element : component) { + uint64_t componentState = internal::getComponentElementState(element); + for (uint64_t choice = choiceIndices[componentState]; choice < choiceIndices[componentState + 1]; ++choice) { + // If the choice is not contained in the component itself, we have to add a similar distribution to the auxiliary state. + if (!internal::componentElementChoicesContains(element, choice)) { rhs.push_back(storm::utility::zero()); - addSspMatrixChoice(choice, _transitionMatrix, inputToSspStateMap, numberOfStatesNotInMecs, currentSspChoice, sspMatrixBuilder); - if (this->isProduceSchedulerSet()) { + addSspMatrixChoice(choice, this->_transitionMatrix, inputToSspStateMap, numberOfNonComponentStates, currentSspChoice, sspMatrixBuilder); + if (sspComponentExitChoicesToOriginalMap) { // Later we need to be able to map this choice back to the original input model - sspMecExitChoicesToOriginalMap.emplace_back(mecState, choice - nondeterministicChoiceIndices[mecState]); + sspComponentExitChoicesToOriginalMap->emplace_back(componentState, choice - choiceIndices[componentState]); } ++currentSspChoice; } } } - // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC. - rhs.push_back(mecLraValues[mecIndex]); - if (this->isProduceSchedulerSet()) { + // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the component. + rhs.push_back(mecLraValues[componentIndex]); + if (sspComponentExitChoicesToOriginalMap) { // Insert some invalid values so we can later detect that this choice is not an exit choice - sspMecExitChoicesToOriginalMap.emplace_back(std::numeric_limits::max(), std::numeric_limits::max()); + sspComponentExitChoicesToOriginalMap->emplace_back(std::numeric_limits::max(), std::numeric_limits::max()); } ++currentSspChoice; } - storm::storage::SparseMatrix sspMatrix = sspMatrixBuilder.build(currentSspChoice, numberOfSspStates, numberOfSspStates); + return std::make_pair(sspMatrixBuilder.build(currentSspChoice, numberOfSspStates, numberOfSspStates), std::move(rhs)); + } + + template + void SparseNondeterministicInfiniteHorizonHelper::constructOptimalChoices(std::vector const& sspChoices, storm::storage::SparseMatrix const& sspMatrix, std::vector const& inputToSspStateMap, storm::storage::BitVector const& statesNotInComponent, uint64_t numberOfNonComponentStates, std::vector> const& sspComponentExitChoicesToOriginalMap) { + // We first take care of non-mec states + storm::utility::vector::setVectorValues(this->_producedOptimalChoices.get(), statesNotInComponent, sspChoices); + // Secondly, we consider MEC states. There are 3 cases for each MEC state: + // 1. The SSP choices encode that we want to stay in the MEC + // 2. The SSP choices encode that we want to leave the MEC and + // a) we take an exit (non-MEC) choice at the given state + // b) we have to take a MEC choice at the given state in a way that eventually an exit state of the MEC is reached + uint64_t exitChoiceOffset = sspMatrix.getRowGroupIndices()[numberOfNonComponentStates]; + for (auto const& mec : *this->_longRunComponentDecomposition) { + // Get the sspState of this MEC (using one representative mec state) + auto const& sspState = inputToSspStateMap[mec.begin()->first]; + uint64_t sspChoiceIndex = sspMatrix.getRowGroupIndices()[sspState] + sspChoices[sspState]; + // Obtain the state and choice of the original model to which the selected choice corresponds. + auto const& originalStateChoice = sspComponentExitChoicesToOriginalMap[sspChoiceIndex - exitChoiceOffset]; + // Check if we are in Case 1 or 2 + if (originalStateChoice.first == std::numeric_limits::max()) { + // The optimal choice is to stay in this mec (Case 1) + // In this case, no further operations are necessary. The scheduler has already been set to the optimal choices during the call of computeLraForMec. + STORM_LOG_ASSERT(sspMatrix.getRow(sspState, sspChoices[sspState]).getNumberOfEntries() == 0, "Expected empty row at choice that stays in MEC."); + } else { + // The best choice is to leave this MEC via the selected state and choice. (Case 2) + // Set the exit choice (Case 2.a) + this->_producedOptimalChoices.get()[originalStateChoice.first] = originalStateChoice.second; + // The remaining states in this MEC need to reach the state with the exit choice with probability 1. (Case 2.b) + // Perform a backwards search from the exit state, only using MEC choices + // We start by setting an invalid choice to all remaining mec states (so that we can easily detect them as unprocessed) + for (auto const& stateActions : mec) { + if (stateActions.first != originalStateChoice.first) { + this->_producedOptimalChoices.get()[stateActions.first] = std::numeric_limits::max(); + } + } + // Ensure that backwards transitions are available + if (this->_backwardTransitions == nullptr) { + this->_computedBackwardTransitions = std::make_unique>(this->_transitionMatrix.transpose(true)); + this->_backwardTransitions = this->_computedBackwardTransitions.get(); + } + // Now start a backwards DFS + std::vector stack = {originalStateChoice.first}; + while (!stack.empty()) { + uint64_t currentState = stack.back(); + stack.pop_back(); + for (auto const& backwardsTransition : this->_backwardTransitions->getRowGroup(currentState)) { + uint64_t predecessorState = backwardsTransition.getColumn(); + if (mec.containsState(predecessorState)) { + auto& selectedPredChoice = this->_producedOptimalChoices.get()[predecessorState]; + if (selectedPredChoice == std::numeric_limits::max()) { + // We don't already have a choice for this predecessor. + // We now need to check whether there is a *MEC* choice leading to currentState + for (auto const& predChoice : mec.getChoicesForState(predecessorState)) { + for (auto const& forwardTransition : this->_transitionMatrix.getRow(predChoice)) { + if (forwardTransition.getColumn() == currentState && !storm::utility::isZero(forwardTransition.getValue())) { + // Playing this choice (infinitely often) will lead to current state (infinitely often)! + selectedPredChoice = predChoice - this->_transitionMatrix.getRowGroupIndices()[predecessorState]; + stack.push_back(predecessorState); + break; + } + } + if (selectedPredChoice != std::numeric_limits::max()) { + break; + } + } + } + } + } + } + } + } + } + + template + std::vector SparseNondeterministicInfiniteHorizonHelper::buildAndSolveSsp(Environment const& env, std::vector const& componentLraValues) { + STORM_LOG_ASSERT(this->_longRunComponentDecomposition != nullptr, "Decomposition not computed, yet."); + + // For fast transition rewriting, we build a mapping from the input state indices to the state indices of a new transition matrix + // which redirects all transitions leading to a former component state to a new auxiliary state. + // There will be one auxiliary state for each component. These states will be appended to the end of the matrix. + + // First gather the states that are part of a component + // and create a mapping from states that lie in a component to the corresponding component index. + storm::storage::BitVector statesInComponents(this->_transitionMatrix.getRowGroupCount()); + std::vector inputToSspStateMap(this->_transitionMatrix.getRowGroupCount(), std::numeric_limits::max()); + for (uint64_t currentComponentIndex = 0; currentComponentIndex < this->_longRunComponentDecomposition->size(); ++currentComponentIndex) { + for (auto const& element : (*this->_longRunComponentDecomposition)[currentComponentIndex]) { + uint64_t state = internal::getComponentElementState(element); + statesInComponents.set(state); + inputToSspStateMap[state] = currentComponentIndex; + } + } + // Now take care of the non-component states. Note that the order of these states will be preserved. + uint64_t numberOfNonComponentStates = 0; + storm::storage::BitVector statesNotInComponent = ~statesInComponents; + for (auto const& nonComponentState : statesNotInComponent) { + inputToSspStateMap[nonComponentState] = numberOfNonComponentStates; + ++numberOfNonComponentStates; + } + // Finalize the mapping for the component states which now still assigns component states to to their component index. + // To make sure that they point to the auxiliary states (located at the end of the SspMatrix), we need to shift them by the + // number of states that are not in a component. + for (auto const& mecState : statesInComponents) { + inputToSspStateMap[mecState] += numberOfNonComponentStates; + } + + // For scheduler extraction, we will need to create a mapping between choices at the auxiliary states and the + // corresponding choices in the original model. + std::vector> sspComponentExitChoicesToOriginalMap; + + // The next step is to create the SSP matrix and the right-hand side of the SSP. + auto sspMatrixVector = buildSspMatrixVector(componentLraValues, inputToSspStateMap, statesNotInComponent, numberOfNonComponentStates, this->isProduceSchedulerSet() ? &sspComponentExitChoicesToOriginalMap : nullptr); // Set-up a solver storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxLinearEquationSolverFactory; storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(env, true, true, this->getOptimizationDirection(), false, this->isProduceSchedulerSet()); requirements.clearBounds(); STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); - std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(env, sspMatrix); + std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(env, sspMatrixVector.first); solver->setHasUniqueSolution(); solver->setHasNoEndComponents(); solver->setTrackScheduler(this->isProduceSchedulerSet()); - auto lowerUpperBounds = std::minmax_element(mecLraValues.begin(), mecLraValues.end()); + auto lowerUpperBounds = std::minmax_element(componentLraValues.begin(), componentLraValues.end()); solver->setLowerBound(*lowerUpperBounds.first); solver->setUpperBound(*lowerUpperBounds.second); solver->setRequirementsChecked(); // Solve the equation system - std::vector x(numberOfSspStates); - solver->solveEquations(env, this->getOptimizationDirection(), x, rhs); + std::vector x(sspMatrixVector.first.getRowGroupCount()); + solver->solveEquations(env, this->getOptimizationDirection(), x, sspMatrixVector.second); // Prepare scheduler (if requested) if (this->isProduceSchedulerSet() && solver->hasScheduler()) { // Translate result for ssp matrix to original model - auto const& sspChoices = solver->getSchedulerChoices(); - // We first take care of non-mec states - storm::utility::vector::setVectorValues(_producedOptimalChoices.get(), statesNotContainedInAnyMec, sspChoices); - // Secondly, we consider MEC states. There are 3 cases for each MEC state: - // 1. The SSP choices encode that we want to stay in the MEC - // 2. The SSP choices encode that we want to leave the MEC and - // a) we take an exit (non-MEC) choice at the given state - // b) we have to take a MEC choice at the given state in a way that eventually an exit state of the MEC is reached - uint64_t exitChoiceOffset = sspMatrix.getRowGroupIndices()[numberOfStatesNotInMecs]; - for (auto const& mec : *_longRunComponentDecomposition) { - // Get the sspState of this MEC (using one representative mec state) - auto const& sspState = inputToSspStateMap[mec.begin()->first]; - uint64_t sspChoiceIndex = sspMatrix.getRowGroupIndices()[sspState] + sspChoices[sspState]; - // Obtain the state and choice of the original model to which the selected choice corresponds. - auto const& originalStateChoice = sspMecExitChoicesToOriginalMap[sspChoiceIndex - exitChoiceOffset]; - // Check if we are in Case 1 or 2 - if (originalStateChoice.first == std::numeric_limits::max()) { - // The optimal choice is to stay in this mec (Case 1) - // In this case, no further operations are necessary. The scheduler has already been set to the optimal choices during the call of computeLraForMec. - STORM_LOG_ASSERT(sspMatrix.getRow(sspState, sspChoices[sspState]).getNumberOfEntries() == 0, "Expected empty row at choice that stays in MEC."); - } else { - // The best choice is to leave this MEC via the selected state and choice. (Case 2) - // Set the exit choice (Case 2.a) - _producedOptimalChoices.get()[originalStateChoice.first] = originalStateChoice.second; - // The remaining states in this MEC need to reach the state with the exit choice with probability 1. (Case 2.b) - // Perform a backwards search from the exit state, only using MEC choices - // We start by setting an invalid choice to all remaining mec states (so that we can easily detect them as unprocessed) - for (auto const& stateActions : mec) { - if (stateActions.first != originalStateChoice.first) { - _producedOptimalChoices.get()[stateActions.first] = std::numeric_limits::max(); - } - } - // Ensure that backwards transitions are available - if (_backwardTransitions == nullptr) { - _computedBackwardTransitions = _transitionMatrix.transpose(true); - _backwardTransitions = &_computedBackwardTransitions; - } - // Now start a backwards DFS - std::vector stack = {originalStateChoice.first}; - while (!stack.empty()) { - uint64_t currentState = stack.back(); - stack.pop_back(); - for (auto const& backwardsTransition : _backwardTransitions->getRowGroup(currentState)) { - uint64_t predecessorState = backwardsTransition.getColumn(); - if (mec.containsState(predecessorState)) { - auto& selectedPredChoice = _producedOptimalChoices.get()[predecessorState]; - if (selectedPredChoice == std::numeric_limits::max()) { - // We don't already have a choice for this predecessor. - // We now need to check whether there is a *MEC* choice leading to currentState - for (auto const& predChoice : mec.getChoicesForState(predecessorState)) { - for (auto const& forwardTransition : _transitionMatrix.getRow(predChoice)) { - if (forwardTransition.getColumn() == currentState && !storm::utility::isZero(forwardTransition.getValue())) { - // Playing this choice (infinitely often) will lead to current state (infinitely often)! - selectedPredChoice = predChoice - nondeterministicChoiceIndices[predecessorState]; - stack.push_back(predecessorState); - break; - } - } - if (selectedPredChoice != std::numeric_limits::max()) { - break; - } - } - } - } - } - } - } - } + constructOptimalChoices(solver->getSchedulerChoices(), sspMatrixVector.first, inputToSspStateMap, statesNotInComponent, numberOfNonComponentStates, sspComponentExitChoicesToOriginalMap); } else { STORM_LOG_ERROR_COND(!this->isProduceSchedulerSet(), "Requested to produce a scheduler, but no scheduler was generated."); } // Prepare result vector. // For efficiency reasons, we re-use the memory of our rhs for this! - std::vector result = std::move(rhs); - result.resize(numberOfStates); + std::vector result = std::move(sspMatrixVector.second); + result.resize(this->_transitionMatrix.getRowGroupCount()); result.shrink_to_fit(); storm::utility::vector::selectVectorValues(result, inputToSspStateMap, x); return result; } - template class SparseNondeterministicInfiniteHorizonHelper; - template class SparseNondeterministicInfiniteHorizonHelper; + template class SparseNondeterministicInfiniteHorizonHelper; + template class SparseNondeterministicInfiniteHorizonHelper; - //template class SparseNondeterministicInfiniteHorizonHelper; - //template class SparseNondeterministicInfiniteHorizonHelper; } } } \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.h b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.h index c11a87a90..5286dfdfc 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.h +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.h @@ -1,20 +1,11 @@ #pragma once -#include "storm/modelchecker/helper/SingleValueModelCheckerHelper.h" +#include "storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.h" namespace storm { - class Environment; - namespace models { - namespace sparse { - template class StandardRewardModel; - } - } namespace storage { - template class Decomposition; - class MaximalEndComponent; - template class SparseMatrix; - class StronglyConnectedComponent; + template class Scheduler; } namespace modelchecker { @@ -23,75 +14,26 @@ namespace storm { /*! * Helper class for model checking queries that depend on the long run behavior of the (nondeterministic) system. * @tparam ValueType the type a value can have - * @tparam Nondeterministic true if there is nondeterminism in the Model (MDP or MA) */ - template - class SparseNondeterministicInfiniteHorizonHelper : public SingleValueModelCheckerHelper { + template + class SparseNondeterministicInfiniteHorizonHelper : public SparseInfiniteHorizonHelper { public: - - /*! - * The type of a component in which the system resides in the long run (BSCC for deterministic models, MEC for nondeterministic models) - */ - using LongRunComponentType = typename std::conditional::type; - /*! * Function mapping from indices to values */ - typedef std::function ValueGetter; + typedef typename SparseInfiniteHorizonHelper::ValueGetter ValueGetter; /*! - * Initializes the helper for a discrete time (i.e. MDP) + * Initializes the helper for a discrete time model (i.e. MDP) */ SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix); /*! - * Initializes the helper for a continuous time (i.e. MA) + * Initializes the helper for a continuous time model (i.e. MA) */ SparseNondeterministicInfiniteHorizonHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates); - /*! - * Provides backward transitions that can be used during the computation. - * Providing them is optional. If they are not provided, they will be computed internally - * Be aware that this class does not take ownership, i.e. the caller has to make sure that the reference to the backwardstransitions remains valid. - */ - void provideBackwardTransitions(storm::storage::SparseMatrix const& backwardsTransitions); - - /*! - * Provides the decomposition into long run components (BSCCs/MECs) that can be used during the computation. - * Providing the decomposition is optional. If it is not provided, they will be computed internally. - * Be aware that this class does not take ownership, i.e. the caller has to make sure that the reference to the decomposition remains valid. - */ - void provideLongRunComponentDecomposition(storm::storage::Decomposition const& decomposition); - - /*! - * Computes the long run average probabilities, i.e., the fraction of the time we are in a psiState - * @return a value for each state - */ - std::vector computeLongRunAverageProbabilities(Environment const& env, storm::storage::BitVector const& psiStates); - - /*! - * Computes the long run average rewards, i.e., the average reward collected per time unit - * @return a value for each state - */ - std::vector computeLongRunAverageRewards(Environment const& env, storm::models::sparse::StandardRewardModel const& rewardModel); - - /*! - * Computes the long run average value given the provided state and action-based rewards. - * @param stateValues a vector containing a value for every state - * @param actionValues a vector containing a value for every choice - * @return a value for each state - */ - std::vector computeLongRunAverageValues(Environment const& env, std::vector const* stateValues = nullptr, std::vector const* actionValues = nullptr); - - /*! - * Computes the long run average value given the provided state and action based rewards - * @param stateValuesGetter a function returning a value for a given state index - * @param actionValuesGetter a function returning a value for a given (global) choice index - * @return a value for each state - */ - std::vector computeLongRunAverageValues(Environment const& env, ValueGetter const& stateValuesGetter, sValueGetter const& actionValuesGetter); - /*! * @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. @@ -116,49 +58,41 @@ namespace storm { * @return the (unique) optimal LRA value for the given component. * @post if scheduler production is enabled and Nondeterministic is true, getProducedOptimalChoices() contains choices for the states of the given component which yield the returned LRA value. Choices for states outside of the component are not affected. */ - template < typename = typename std::enable_if< true >::type > - ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& component); - template < typename = typename std::enable_if< false >::type > - ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& component); + virtual ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::MaximalEndComponent const& component) override; protected: - /*! - * @return true iff this is a computation on a continuous time model (i.e. MA) - */ - bool isContinuousTime() const; + virtual void createDecomposition() override; - /*! - * Checks if the component can trivially be solved without much overhead. - * @return either true and the (unique) optimal LRA value for the given component or false and an arbitrary value - * @post if scheduler production is enabled and Nondeterministic is true, getProducedOptimalChoices() contains choices for the states of the given component which yield the returned LRA value. Choices for states outside of the component are not affected. - */ - std::pair computeLraForTrivialComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& component); + std::pair computeLraForTrivialMec(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::MaximalEndComponent const& mec); /*! * As computeLraForMec but uses value iteration as a solution method (independent of what is set in env) */ - ValueType computeLraForMecVi(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& mec); + ValueType computeLraForMecVi(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::MaximalEndComponent const& mec); + /*! * As computeLraForMec but uses linear programming as a solution method (independent of what is set in env) * @see Guck et al.: Modelling and Analysis of Markov Reward Automata (ATVA'14), https://doi.org/10.1007/978-3-319-11936-6_13 */ - ValueType computeLraForMecLp(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, LongRunComponentType const& mec); + ValueType computeLraForMecLp(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::MaximalEndComponent const& mec); + + std::pair, std::vector> buildSspMatrixVector(std::vector const& mecLraValues, std::vector const& inputToSspStateMap, storm::storage::BitVector const& statesNotInComponent, uint64_t numberOfNonComponentStates, std::vector>* sspComponentExitChoicesToOriginalMap); + + /*! + * @pre a choice for each state within a component must be set such that the choices yield optimal values w.r.t. the individual components. + * Translates optimal choices for MECS and SSP to the original model. + * @post getProducedOptimalChoices() contains choices for all input model states which yield the returned LRA values. + */ + void constructOptimalChoices(std::vector const& sspChoices, storm::storage::SparseMatrix const& sspMatrix, std::vector const& inputToSspStateMap, storm::storage::BitVector const& statesNotInComponent, uint64_t numberOfNonComponentStates, std::vector> const& sspComponentExitChoicesToOriginalMap); /*! + * @pre if scheduler production is enabled a choice for each state within a component must be set such that the choices yield optimal values w.r.t. the individual components. * @return Lra values for each state + * @post if scheduler production is enabled getProducedOptimalChoices() contains choices for all input model states which yield the returned LRA values. */ - std::vector buildAndSolveSsp(Environment const& env, std::vector const& mecLraValues); + virtual std::vector buildAndSolveSsp(Environment const& env, std::vector const& mecLraValues) override; - private: - storm::storage::SparseMatrix const& _transitionMatrix; - storm::storage::SparseMatrix const* _backwardTransitions; - std::unique_ptr> _computedBackwardTransitions; - storm::storage::Decomposition const* _longRunComponentDecomposition; - std::unique_ptr> _computedLongRunComponentDecomposition; - storm::storage::BitVector const* _markovianStates; - std::vector const* _exitRates; - boost::optional> _producedOptimalChoices; };