From 6e55dba8d4d11491ccf018bcc317c3f7517e361a Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Wed, 5 Aug 2020 11:39:03 +0200 Subject: [PATCH] Moved LraViHelper to a separate file. Merged MDP and MA implementation. --- .../infinitehorizon/internal/LraViHelper.cpp | 517 ++++++++++++++++++ .../infinitehorizon/internal/LraViHelper.h | 151 +++++ 2 files changed, 668 insertions(+) create mode 100644 src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp create mode 100644 src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h diff --git a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp new file mode 100644 index 000000000..d4c3ba856 --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp @@ -0,0 +1,517 @@ +#include "LraViHelper.h" + +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/solver/Multiplier.h" + +#include "storm/storage/MaximalEndComponent.h" +#include "storm/storage/StronglyConnectedComponent.h" + +#include "storm/utility/graph.h" +#include "storm/utility/vector.h" +#include "storm/utility/macros.h" +#include "storm/utility/SignalHandler.h" + +#include "storm/environment/solver/SolverEnvironment.h" +#include "storm/environment/solver/LongRunAverageSolverEnvironment.h" +#include "storm/environment/solver/MinMaxSolverEnvironment.h" + + +#include "storm/exceptions/UnmetRequirementException.h" + + +namespace storm { + namespace modelchecker { + namespace helper { + namespace internal { + + /// Auxiliary functions that deal with the different kinds of components (MECs on potentially nondeterministic models and BSCCs on deterministic models) + // BSCCS: + uint64_t getComponentElementState(typename storm::storage::StronglyConnectedComponent::value_type const& element) { return element; } + uint64_t getComponentElementChoiceCount(typename storm::storage::StronglyConnectedComponent::value_type const& element) { return 1; } // Assumes deterministic model! + uint64_t const* getComponentChoicesBegin(typename storm::storage::StronglyConnectedComponent::value_type const& element) { return &element; } + uint64_t const* getComponentChoicesEnd(typename storm::storage::StronglyConnectedComponent::value_type const& element) { return &element + 1; } + // MECS: + uint64_t getComponentElementState(typename storm::storage::MaximalEndComponent::map_type::value_type const& element) { return element.first; } + uint64_t getComponentElementChoiceCount(typename storm::storage::MaximalEndComponent::map_type::value_type const& element) { return element.second.size(); } + typename storm::storage::MaximalEndComponent::set_type::const_iterator getComponentChoicesBegin(typename storm::storage::MaximalEndComponent::map_type::value_type const& element) { return element.second.begin(); } + typename storm::storage::MaximalEndComponent::set_type::const_iterator getComponentChoicesEnd(typename storm::storage::MaximalEndComponent::map_type::value_type const& element) { return element.second.end(); } + + template + LraViHelper::LraViHelper(ComponentType const& component, storm::storage::SparseMatrix const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates, std::vector const* exitRates) : _component(component), _transitionMatrix(transitionMatrix), _timedStates(timedStates), _hasInstantStates(TransitionsType == LraViTransitionsType::DetTsNondetIs || TransitionsType == LraViTransitionsType::DetTsDetIs) { + // Run through the component and collect some data: + // We create two submodels, one consisting of the timed states of the component and one consisting of the instant states of the component. + // For this, we create a state index map that point from state indices of the input model to indices of the corresponding submodel of that state. + boost::container::flat_map toSubModelStateMapping; + // We also obtain state and choices counts of the two submodels + uint64_t numTsSubModelStates(0), numTsSubModelChoices(0); + uint64_t numIsSubModelStates(0), numIsSubModelChoices(0); + // We will need to uniformize the timed MEC states by introducing a selfloop. + // For this, we need to find a uniformization rate which will be a little higher (given by aperiodicFactor) than the maximum rate occurring in the component. + _uniformizationRate = exitRates == nullptr ? storm::utility::one() : storm::utility::zero(); + // Now run over the MEC and collect the required data. + for (auto const& element : _component) { + uint64_t const& componentState = getComponentElementState(element); + if (isTimedState(componentState)) { + toSubModelStateMapping.emplace(componentState, numTsSubModelStates); + ++numTsSubModelStates; + numTsSubModelChoices += getComponentElementChoiceCount(element); + STORM_LOG_ASSERT(nondetTs() || getComponentElementChoiceCount(element) == 1, "Timed state has multiple choices but only a single choice was expected."); + if (exitRates) { + _uniformizationRate = std::max(_uniformizationRate, (*exitRates)[componentState]); + } + } else { + toSubModelStateMapping.emplace(componentState, numIsSubModelStates); + ++numIsSubModelStates; + numIsSubModelChoices += getComponentElementChoiceCount(element); + STORM_LOG_ASSERT(nondetIs() || getComponentElementChoiceCount(element) == 1, "Instant state has multiple choices but only a single choice was expected."); + } + } + assert(numIsSubModelStates + numTsSubModelStates == _component.size()); + assert(_hasInstantStates || numIsSubModelStates == 0); + STORM_LOG_ASSERT(nondetTs() || numTsSubModelStates == numTsSubModelChoices, "Unexpected choice count of deterministic timed submodel."); + STORM_LOG_ASSERT(nondetIs() || numIsSubModelStates == numIsSubModelChoices, "Unexpected choice count of deterministic instant submodel."); + _hasInstantStates = _hasInstantStates && numIsSubModelStates > 0; + STORM_LOG_THROW(numTsSubModelStates > 0, storm::exceptions::InvalidOperationException, "Bottom Component has no timed states. Computation of Long Run Average values not supported. Is this a Markov Automaton with Zeno behavior?"); + + // We make sure that every timed state gets a selfloop to make the model aperiodic + _uniformizationRate *= storm::utility::one() + aperiodicFactor; + + // Now build the timed and the instant submodels. + // In addition, we also need the transitions between the two. + storm::storage::SparseMatrixBuilder tsTransitionsBuilder(numTsSubModelChoices, numTsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0); + storm::storage::SparseMatrixBuilder tsToIsTransitionsBuilder, isTransitionsBuilder, isToTsTransitionsBuilder; + if (_hasInstantStates) { + tsToIsTransitionsBuilder = storm::storage::SparseMatrixBuilder(numTsSubModelChoices, numIsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0); + isTransitionsBuilder = storm::storage::SparseMatrixBuilder(numIsSubModelChoices, numIsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0); + isToTsTransitionsBuilder = storm::storage::SparseMatrixBuilder(numIsSubModelChoices, numTsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0); + _IsChoiceValues.reserve(numIsSubModelChoices); + } + ValueType uniformizationFactor = storm::utility::one() / _uniformizationRate; + uint64_t currTsRow = 0; + uint64_t currIsRow = 0; + for (auto const& element : _component) { + uint64_t const& componentState = getComponentElementState(element); + if (isTimedState(componentState)) { + // The currently processed state is timed. + if (nondetTs()) { + tsTransitionsBuilder.newRowGroup(currTsRow); + tsToIsTransitionsBuilder.newRowGroup(currTsRow); + } + // We need to uniformize which means that a diagonal entry for the selfloop will be inserted. + // If there are exit rates, the uniformization factor needs to be updated. + if (exitRates) { + uniformizationFactor = (*exitRates)[componentState] / _uniformizationRate; + } + ValueType selfLoopProb = storm::utility::one() - uniformizationFactor; + uint64_t selfLoopColumn = toSubModelStateMapping[componentState]; + for (auto componentChoiceIt = getComponentChoicesBegin(element); componentChoiceIt != getComponentChoicesEnd(element); ++componentChoiceIt) { + bool insertedDiagElement = false; + for (auto const& entry : this->_transitionMatrix.getRow(*componentChoiceIt)) { + uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()]; + if (isTimedState(entry.getColumn())) { + // We have a transition from a timed state to a timed state + STORM_LOG_ASSERT(subModelColumn < numTsSubModelStates, "Invalid state for timed submodel"); + if (!insertedDiagElement && subModelColumn > selfLoopColumn) { + // We passed the diagonal entry, so add it now before moving on to the next entry + tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb); + insertedDiagElement = true; + } + if (!insertedDiagElement && subModelColumn == selfLoopColumn) { + // The current entry is the diagonal (selfloop) entry + tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb + uniformizationFactor * entry.getValue()); + insertedDiagElement = true; + } else { + // The diagonal element either has been inserted already or still lies in front + tsTransitionsBuilder.addNextValue(currTsRow, subModelColumn, uniformizationFactor * entry.getValue()); + } + } else { + // We have a transition from a timed to a instant state + STORM_LOG_ASSERT(subModelColumn < numIsSubModelStates, "Invalid state for instant submodel"); + tsToIsTransitionsBuilder.addNextValue(currTsRow, subModelColumn, uniformizationFactor * entry.getValue()); + } + } + // If the diagonal entry for the MS matrix still has not been set, we do that now + if (!insertedDiagElement) { + tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb); + } + ++currTsRow; + } + } else { + // The currently processed state is instant + if (nondetIs()) { + isTransitionsBuilder.newRowGroup(currIsRow); + isToTsTransitionsBuilder.newRowGroup(currIsRow); + } + for (auto componentChoiceIt = getComponentChoicesBegin(element); componentChoiceIt != getComponentChoicesEnd(element); ++componentChoiceIt) { + for (auto const& entry : this->_transitionMatrix.getRow(*componentChoiceIt)) { + uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()]; + if (isTimedState(entry.getColumn())) { + // We have a transition from an instant state to a timed state + STORM_LOG_ASSERT(subModelColumn < numTsSubModelStates, "Invalid state for timed submodel"); + isToTsTransitionsBuilder.addNextValue(currIsRow, subModelColumn, entry.getValue()); + } else { + // We have a transition from an instant to an instant state + STORM_LOG_ASSERT(subModelColumn < numIsSubModelStates, "Invalid state for instant submodel"); + isTransitionsBuilder.addNextValue(currIsRow, subModelColumn, entry.getValue()); + } + } + ++currIsRow; + } + } + } + _TsTransitions = tsTransitionsBuilder.build(); + if (_hasInstantStates) { + _TsToIsTransitions = tsToIsTransitionsBuilder.build(); + _IsTransitions = isTransitionsBuilder.build(); + _IsToTsTransitions = isToTsTransitionsBuilder.build(); + } + } + + template + ValueType LraViHelper::performValueIteration(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates, storm::solver::OptimizationDirection const* dir, std::vector* choices) { + + initializeNewValues(stateValueGetter, actionValueGetter, exitRates); + + ValueType precision = storm::utility::convertNumber(env.solver().lra().getPrecision()); + bool relative = env.solver().lra().getRelativeTerminationCriterion(); + boost::optional maxIter; + if (env.solver().lra().isMaximalIterationCountSet()) { + maxIter = env.solver().lra().getMaximalIterationCount(); + } + + // start the iterations + ValueType result = storm::utility::zero(); + uint64_t iter = 0; + while (!maxIter.is_initialized() || iter < maxIter.get()) { + ++iter; + performIterationStep(env, dir); + + // Check if we are done + auto convergenceCheckResult = checkConvergence(relative, precision); + result = convergenceCheckResult.currentValue; + if (convergenceCheckResult.isPrecisionAchieved) { + break; + } + if (storm::utility::resources::isTerminate()) { + break; + } + // If there will be a next iteration, we have to prepare it. + prepareNextIteration(env); + + } + if (maxIter.is_initialized() && iter == maxIter.get()) { + STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations."); + } else if (storm::utility::resources::isTerminate()) { + STORM_LOG_WARN("LRA computation aborted after " << iter << " iterations."); + } else { + STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations."); + } + + if (choices) { + // We will be doing one more iteration step and track scheduler choices this time. + prepareNextIteration(env); + performIterationStep(env, dir, choices); + } + return result; + } + + template + void LraViHelper::initializeNewValues(ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates) { + // clear potential old values and reserve enough space for new values + _TsChoiceValues.clear(); + _TsChoiceValues.reserve(_TsTransitions.getRowCount()); + if (_hasInstantStates) { + _IsChoiceValues.clear(); + _IsChoiceValues.reserve(_IsTransitions.getRowCount()); + } + + // Set the new choice-based values + ValueType actionRewardScalingFactor = storm::utility::one() / _uniformizationRate; + for (auto const& element : _component) { + uint64_t const& componentState = getComponentElementState(element); + if (isTimedState(componentState)) { + if (exitRates) { + actionRewardScalingFactor = (*exitRates)[componentState] / _uniformizationRate; + } + for (auto componentChoiceIt = getComponentChoicesBegin(element); componentChoiceIt != getComponentChoicesEnd(element); ++componentChoiceIt) { + // Compute the values obtained for this choice. + _TsChoiceValues.push_back(stateValueGetter(componentState) / _uniformizationRate + actionValueGetter(*componentChoiceIt) * actionRewardScalingFactor); + } + } else { + for (auto componentChoiceIt = getComponentChoicesBegin(element); componentChoiceIt != getComponentChoicesEnd(element); ++componentChoiceIt) { + // Compute the values obtained for this choice. + // State values do not count here since no time passes in instant states. + _IsChoiceValues.push_back(actionValueGetter(*componentChoiceIt)); + } + } + } + + // Set-up new iteration vectors for timed states + _Tsx1.assign(_TsTransitions.getRowGroupCount(), storm::utility::zero()); + _Tsx2 = _Tsx1; + + if (_hasInstantStates) { + // Set-up vectors for storing intermediate results for instant states. + _Isx.resize(_IsTransitions.getRowGroupCount(), storm::utility::zero()); + _Isb = _IsChoiceValues; + } + } + + template + void LraViHelper::prepareSolversAndMultipliers(const Environment& env, storm::solver::OptimizationDirection const* dir) { + _TsMultiplier = storm::solver::MultiplierFactory().create(env, _TsTransitions); + if (_hasInstantStates) { + if (_IsTransitions.getNonzeroEntryCount() > 0) { + // Set-up a solver for transitions within instant states + _IsSolverEnv = std::make_unique(env); + if (env.solver().isForceSoundness()) { + // To get correct results, the inner equation systems are solved exactly. + // TODO investigate how an error would propagate + _IsSolverEnv->solver().setForceExact(true); + } + bool isAcyclic = !storm::utility::graph::hasCycle(_IsTransitions); + if (isAcyclic) { + STORM_LOG_INFO("Instant transitions are acyclic."); + if (_IsSolverEnv->solver().minMax().isMethodSetFromDefault()) { + _IsSolverEnv->solver().minMax().setMethod(storm::solver::MinMaxMethod::Acyclic); + } + if (_IsSolverEnv->solver().isLinearEquationSolverTypeSetFromDefaultValue()) { + _IsSolverEnv->solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Acyclic); + } + } + if (nondetIs()) { + storm::solver::GeneralMinMaxLinearEquationSolverFactory factory; + _NondetIsSolver = factory.create(*_IsSolverEnv, _IsTransitions); + _NondetIsSolver->setHasUniqueSolution(true); // Assume non-zeno MA + _NondetIsSolver->setHasNoEndComponents(true); // assume non-zeno MA + _NondetIsSolver->setCachingEnabled(true); + auto req = _NondetIsSolver->getRequirements(*_IsSolverEnv, *dir); + req.clearUniqueSolution(); + if (isAcyclic) { + req.clearAcyclic(); + } + // Computing a priori lower/upper bounds is not particularly easy, as there might be selfloops with high probabilities + // Which accumulate a lot of reward. Moreover, the right-hand-side of the equation system changes dynamically. + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "The solver requirement " << req.getEnabledRequirementsAsString() << " has not been cleared."); + _NondetIsSolver->setRequirementsChecked(true); + } else { + storm::solver::GeneralLinearEquationSolverFactory factory; + if (factory.getEquationProblemFormat(*_IsSolverEnv) != storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem) { + // We need to convert the transition matrix connecting instant states + // TODO: This could have been done already during construction of the matrix. + // Insert diagonal entries. + storm::storage::SparseMatrix converted(_IsTransitions, true); + // Compute A' = 1-A + converted.convertToEquationSystem(); + STORM_LOG_WARN("The selected equation solver requires to create a temporary " << converted.getDimensionsAsString()); + // Note that the solver has ownership of the converted matrix. + _DetIsSolver = factory.create(*_IsSolverEnv, std::move(converted)); + } else { + _DetIsSolver = factory.create(*_IsSolverEnv, _IsTransitions); + } + _DetIsSolver->setCachingEnabled(true); + auto req = _DetIsSolver->getRequirements(*_IsSolverEnv); + if (isAcyclic) { + req.clearAcyclic(); + } + // A priori lower/upper bounds are hard (see MinMax version of this above) + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "The solver requirement " << req.getEnabledRequirementsAsString() << " has not been cleared."); + } + } + + // Set up multipliers for transitions connecting timed and instant states + _TsToIsMultiplier = storm::solver::MultiplierFactory().create(env, _TsToIsTransitions); + _IsToTsMultiplier = storm::solver::MultiplierFactory().create(env, _IsToTsTransitions); + } + } + + template + void LraViHelper::setInputModelChoices(std::vector& choices, std::vector const& localMecChoices, bool setChoiceZeroToTimedStates, bool setChoiceZeroToInstantStates) const { + // Transform the local choices (within this mec) to choice indices for the input model + uint64_t localState = 0; + for (auto const& element : _component) { + uint64_t elementState = getComponentElementState(element); + if ((setChoiceZeroToTimedStates && isTimedState(elementState)) || (setChoiceZeroToInstantStates && !isTimedState(elementState))) { + choices[elementState] = 0; + } else { + uint64_t choice = localMecChoices[localState]; + STORM_LOG_ASSERT(choice < getComponentElementChoiceCount(element), "The selected choice does not seem to exist."); + uint64_t globalChoiceIndex = *(getComponentChoicesBegin(element) + choice); + choices[elementState] = globalChoiceIndex - _transitionMatrix.getRowGroupIndices()[elementState]; + ++localState; + } + } + STORM_LOG_ASSERT(localState == localMecChoices.size(), "Did not traverse all component states."); + } + + template + void LraViHelper::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const* dir, std::vector* choices) { + STORM_LOG_ASSERT(!((nondetTs() || nondetIs()) && dir == nullptr), "No optimization direction provided for model with nondeterminism"); + // Initialize value vectors, multiplers, and solver if this has not been done, yet + if (!_TsMultiplier) { + prepareSolversAndMultipliers(env, dir); + } + + // Compute new x values for the timed states + // Flip what is new and what is old + _Tsx1IsCurrent = !_Tsx1IsCurrent; + // At this point, xOld() points to what has been computed in the most recent call of performIterationStep (initially, this is the 0-vector). + // The result of this ongoing computation will be stored in xNew() + + // Compute the values obtained by a single uniformization step between timed states only + if (nondetTs()) { + if (choices == nullptr) { + _TsMultiplier->multiplyAndReduce(env, *dir, xOld(), &_TsChoiceValues, xNew()); + } else { + // Also keep track of the choices made. + std::vector tsChoices(_TsTransitions.getRowGroupCount()); + _TsMultiplier->multiplyAndReduce(env, *dir, xOld(), &_TsChoiceValues, xNew(), &tsChoices); + // Note that nondeterminism within the timed states means that there can not be instant states (We either have MDPs or MAs) + // Hence, in this branch we don't have to care for choices at instant states. + STORM_LOG_ASSERT(!_hasInstantStates, "Nondeterministic timed states are only supported if there are no instant states."); + setInputModelChoices(*choices, tsChoices); + } + } else { + _TsMultiplier->multiply(env, xOld(), &_TsChoiceValues, xNew()); + } + if (_hasInstantStates) { + // Add the values obtained by taking a single uniformization step that leads to an instant state followed by arbitrarily many instant steps. + // First compute the total values when taking arbitrarily many instant transitions (in no time) + if (_NondetIsSolver) { + // We might need to track the optimal choices. + if (choices == nullptr) { + _NondetIsSolver->solveEquations(*_IsSolverEnv, *dir, _Isx, _Isb); + } else { + _NondetIsSolver->setTrackScheduler(); + _NondetIsSolver->solveEquations(*_IsSolverEnv, *dir, _Isx, _Isb); + setInputModelChoices(*choices, _NondetIsSolver->getSchedulerChoices(), true); + } + } else if (_DetIsSolver) { + _DetIsSolver->solveEquations(*_IsSolverEnv, _Isx, _Isb); + } else { + STORM_LOG_ASSERT(_IsTransitions.getNonzeroEntryCount() == 0, "If no solver was initialized, an empty matrix would have been expected."); + if (nondetIs()) { + if (choices == nullptr) { + storm::utility::vector::reduceVectorMinOrMax(*dir, _Isb, _Isx, _IsTransitions.getRowGroupIndices()); + } else { + std::vector psChoices(_IsTransitions.getRowGroupCount()); + storm::utility::vector::reduceVectorMinOrMax(*dir, _Isb, _Isx, _IsTransitions.getRowGroupIndices(), &psChoices); + setInputModelChoices(*choices, psChoices, true); + } + } else { + // For deterministic instant states, there is nothing to reduce, i.e., we could just set _Isx = _Isb. + // For efficiency reasons, we do a swap instead: + _Isx.swap(_Isb); + // Note that at this point we have changed the contents of _Isb, but they will be overwritten anyway. + if (choices) { + // Set choice 0 to all states. + setInputModelChoices(*choices, {}, true, true); + } + } + } + // Now add the (weighted) values of the instant states to the values of the timed states. + _TsToIsMultiplier->multiply(env, _Isx, &xNew(), xNew()); + } + } + + template + typename LraViHelper::ConvergenceCheckResult LraViHelper::checkConvergence(bool relative, ValueType precision) const { + STORM_LOG_ASSERT(_TsMultiplier, "tried to check for convergence without doing an iteration first."); + // All values are scaled according to the uniformizationRate. + // We need to 'revert' this scaling when computing the absolute precision. + // However, for relative precision, the scaling cancels out. + ValueType threshold = relative ? precision : ValueType(precision / _uniformizationRate); + + ConvergenceCheckResult res = { true, storm::utility::one() }; + // Now check whether the currently produced results are precise enough + STORM_LOG_ASSERT(threshold > storm::utility::zero(), "Did not expect a non-positive threshold."); + auto x1It = xOld().begin(); + auto x1Ite = xOld().end(); + auto x2It = xNew().begin(); + ValueType maxDiff = (*x2It - *x1It); + ValueType minDiff = maxDiff; + // The difference between maxDiff and minDiff is zero at this point. Thus, it doesn't make sense to check the threshold now. + for (++x1It, ++x2It; x1It != x1Ite; ++x1It, ++x2It) { + ValueType diff = (*x2It - *x1It); + // Potentially update maxDiff or minDiff + bool skipCheck = false; + if (maxDiff < diff) { + maxDiff = diff; + } else if (minDiff > diff) { + minDiff = diff; + } else { + skipCheck = true; + } + // Check convergence + if (!skipCheck && (maxDiff - minDiff) > (relative ? (threshold * minDiff) : threshold)) { + res.isPrecisionAchieved = false; + break; + } + } + + // Compute the average of the maximal and the minimal difference. + ValueType avgDiff = (maxDiff + minDiff) / (storm::utility::convertNumber(2.0)); + + // "Undo" the scaling of the values + res.currentValue = avgDiff * _uniformizationRate; + return res; + } + + template + void LraViHelper::prepareNextIteration(Environment const& env) { + // To avoid large (and numerically unstable) x-values, we substract a reference value. + ValueType referenceValue = xNew().front(); + storm::utility::vector::applyPointwise(xNew(), xNew(), [&referenceValue] (ValueType const& x_i) -> ValueType { return x_i - referenceValue; }); + if (_hasInstantStates) { + // Update the RHS of the equation system for the instant states by taking the new values of timed states into account. + STORM_LOG_ASSERT(!nondetTs(), "Nondeterministic timed states not expected when there are also instant states."); + _IsToTsMultiplier->multiply(env, xNew(), &_IsChoiceValues, _Isb); + } + } + + template + bool LraViHelper::isTimedState(uint64_t const& inputModelStateIndex) const { + STORM_LOG_ASSERT(!_hasInstantStates || inputModelStateIndex < _timedStates->size(), "Unable to determine whether state " << inputModelStateIndex << " is timed."); + return !_hasInstantStates || _timedStates->get(inputModelStateIndex); + } + + template + std::vector& LraViHelper::xNew() { + return _Tsx1IsCurrent ? _Tsx1 : _Tsx2; + } + + template + std::vector const& LraViHelper::xNew() const { + return _Tsx1IsCurrent ? _Tsx1 : _Tsx2; + } + + template + std::vector& LraViHelper::xOld() { + return _Tsx1IsCurrent ? _Tsx2 : _Tsx1; + } + + template + std::vector const& LraViHelper::xOld() const { + return _Tsx1IsCurrent ? _Tsx2 : _Tsx1; + } + + template + bool LraViHelper::nondetTs() const { + return TransitionsType == LraViTransitionsType::NondetTsNoIs; + } + + template + bool LraViHelper::nondetIs() const { + return TransitionsType == LraViTransitionsType::DetTsNondetIs; + } + + template class LraViHelper; + template class LraViHelper; + template class LraViHelper; + template class LraViHelper; + + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h new file mode 100644 index 000000000..3fb5066c4 --- /dev/null +++ b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h @@ -0,0 +1,151 @@ +#pragma once + + +#include "storm/storage/SparseMatrix.h" + +namespace storm { + class Environment; + + namespace solver { + template class LinearEquationSolver; + template class MinMaxLinearEquationSolver; + template class Multiplier; + } + + namespace modelchecker { + namespace helper { + namespace internal { + + /*! + * Specifies differnt kinds of transition types with which this helper can be used + * Ts means timed states (cf. Markovian states in a Markov Automaton) and Is means instant states (cf. probabilistic states in a Markov automaton). + * The way to think about this is that time can only pass in a timed state, whereas transitions emerging from an instant state fire immediately + * In an MDP, all states are seen as timed. + * In this enum, we also specify whether there can be a nondeterministic choice at the corresponding states or not. + */ + enum class LraViTransitionsType { + DetTsNoIs, /// deterministic choice at timed states, no instant states (as in DTMCs and CTMCs) + DetTsNondetIs, /// deterministic choice at timed states, nondeterministic choice at instant states (as in Markov Automata) + DetTsDetIs, /// deterministic choice at timed states, deterministic choice at instant states (as in Markov Automata without any nondeterminisim) + NondetTsNoIs /// nondeterministic choice at timed states, no instant states (as in MDPs) + }; + + /*! + * Helper class that performs iterations of the value iteration method. + * The purpose of the template parameters ComponentType and TransitionsType are used to make this work for various model types. + * + * @see Ashok et al.: Value Iteration for Long-Run Average Reward in Markov Decision Processes (CAV'17), https://doi.org/10.1007/978-3-319-63387-9_10 + * @see Butkova, Wimmer, Hermanns: Long-Run Rewards for Markov Automata (TACAS'17), https://doi.org/10.1007/978-3-662-54580-5_11 + * + * @tparam ValueType The type of a value + * @tparam ComponentType The type of a 'bottom component' of the model (e.g. a BSCC (for purely deterministic models) or a MEC (for models with potential nondeterminism). + * @tparam TransitionsType The kind of transitions that occur. + */ + template + class LraViHelper { + public: + /// Function mapping from indices to values + typedef std::function ValueGetter; + + /*! + * Initializes a new VI helper for the provided MEC or BSCC + * @param component the MEC or BSCC + * @param transitionMatrix The transition matrix of the input model + * @param aperiodicFactor a non-zero factor that is used for making the MEC aperiodic (by adding selfloops to each state) + * @param timedStates States in which time can pass (Markovian states in a Markov automaton). If nullptr, it is assumed that all states are timed states + * @param exitRates The exit rates of the timed states (relevant for continuous time models). If nullptr, all rates are assumed to be 1 (which corresponds to a discrete time model) + * @note All indices and vectors must be w.r.t. the states as described by the provided transition matrix + */ + LraViHelper(ComponentType const& component, storm::storage::SparseMatrix const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates = nullptr, std::vector const* exitRates = nullptr); + + /*! + * Performs value iteration with the given state- and action values. + * @param env The environment, containing information on the precision of this computation. + * @param stateValueGetter function that returns for each state index (w.r.t. the input transition matrix) the reward for staying in state. Will only be called for timed states. + * @param actionValueGetter function that returns for each global choice index (w.r.t. the input transition matrix) the reward for taking that choice + * @param exitRates (as in the constructor) + * @param dir Optimization direction. Must be not nullptr in case of nondeterminism + * @param choices if not nullptr, the optimal choices will be inserted in this vector. The vector's size must then be equal to the number of row groups of the input transition matrix. + * @return The (optimal) long run average value of the specified component. + * @note it is possible to call this method multiple times with different values. However, other changes to the environment or the optimization direction might not have the expected effect due to caching. + */ + ValueType performValueIteration(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates = nullptr, storm::solver::OptimizationDirection const* dir = nullptr, std::vector* choices = nullptr); + + private: + + /*! + * Initializes the value iterations with the provided values. + * Resets all information from potential previous calls. + * Must be called before the first call to performIterationStep. + * @param stateValueGetter Function that returns for each state index (w.r.t. the input transitions) the value (e.g. reward) for that state + * @param stateValueGetter Function that returns for each global choice index (w.r.t. the input transitions) the value (e.g. reward) for that choice + */ + void initializeNewValues(ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates = nullptr); + + /*! + * Performs a single iteration step. + * @param env The environment. + * @param dir The optimization direction. Has to be given if there is nondeterminism (otherwise it will be ignored) + * @param choices If given, the optimal choices will be inserted at the appropriate states. + * Note that these choices will be inserted w.r.t. the original model states/choices, i.e. the size of the vector should match the state-count of the input model + * @pre when calling this the first time, initializeNewValues must have been called before. Moreover, prepareNextIteration must be called between two calls of this. + */ + void performIterationStep(Environment const& env, storm::solver::OptimizationDirection const* dir = nullptr, std::vector* choices = nullptr); + + struct ConvergenceCheckResult { + bool isPrecisionAchieved; + ValueType currentValue; + }; + + /*! + * Checks whether the curently computed value achieves the desired precision + */ + ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) const; + + /*! + * Must be called between two calls of performIterationStep. + */ + void prepareNextIteration(Environment const& env); + + /// Prepares the necessary solvers and multipliers for doing the iterations. + void prepareSolversAndMultipliers(Environment const& env, storm::solver::OptimizationDirection const* dir = nullptr); + + void setInputModelChoices(std::vector& choices, std::vector const& localMecChoices, bool setChoiceZeroToMarkovianStates = false, bool setChoiceZeroToProbabilisticStates = false) const; + + /// Returns true iff the given state is a timed state + bool isTimedState(uint64_t const& inputModelStateIndex) const; + + /// The result for timed states of the most recent iteration + std::vector& xNew(); + std::vector const& xNew() const; + + /// The result for timed states of the previous iteration + std::vector& xOld(); + std::vector const& xOld() const; + + /// @return true iff there potentially is a nondeterministic choice at timed states + bool nondetTs() const; + + /// @return true iff there potentially is a nondeterministic choice at instant states. Returns false if there are no instant states. + bool nondetIs() const; + + + + ComponentType const& _component; + storm::storage::SparseMatrix const& _transitionMatrix; + storm::storage::BitVector const* _timedStates; // e.g. Markovian states of a Markov automaton. + bool _hasInstantStates; + ValueType _uniformizationRate; + storm::storage::SparseMatrix _TsTransitions, _TsToIsTransitions, _IsTransitions, _IsToTsTransitions; + std::vector _Tsx1, _Tsx2, _TsChoiceValues; + bool _Tsx1IsCurrent; + std::vector _Isx, _Isb, _IsChoiceValues; + std::unique_ptr> _TsMultiplier, _TsToIsMultiplier, _IsToTsMultiplier; + std::unique_ptr> _NondetIsSolver; + std::unique_ptr> _DetIsSolver; + std::unique_ptr _IsSolverEnv; + }; + } + } + } +} \ No newline at end of file