sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Moved LraViHelper to a separate file. Merged MDP and MA implementation.

tempestpy_adaptions
Tim Quatmann 4 years ago
parent
a618147192
commit
6e55dba8d4
2 changed files with 668 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 517
      src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp
  2. 151
      src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h

517
src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp
View File

@ -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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
LraViHelper<ValueType, ComponentType, TransitionsType>::LraViHelper(ComponentType const& component, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates, std::vector<ValueType> 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<uint64_t, uint64_t> 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<ValueType>() : storm::utility::zero<ValueType>();
// 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<ValueType>() + aperiodicFactor;
// Now build the timed and the instant submodels.
// In addition, we also need the transitions between the two.
storm::storage::SparseMatrixBuilder<ValueType> tsTransitionsBuilder(numTsSubModelChoices, numTsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0);
storm::storage::SparseMatrixBuilder<ValueType> tsToIsTransitionsBuilder, isTransitionsBuilder, isToTsTransitionsBuilder;
if (_hasInstantStates) {
tsToIsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numTsSubModelChoices, numIsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0);
isTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numIsSubModelChoices, numIsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0);
isToTsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numIsSubModelChoices, numTsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0);
_IsChoiceValues.reserve(numIsSubModelChoices);
}
ValueType uniformizationFactor = storm::utility::one<ValueType>() / _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<ValueType>() - 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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
ValueType LraViHelper<ValueType, ComponentType, TransitionsType>::performValueIteration(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector<ValueType> const* exitRates, storm::solver::OptimizationDirection const* dir, std::vector<uint64_t>* choices) {
initializeNewValues(stateValueGetter, actionValueGetter, exitRates);
ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision());
bool relative = env.solver().lra().getRelativeTerminationCriterion();
boost::optional<uint64_t> maxIter;
if (env.solver().lra().isMaximalIterationCountSet()) {
maxIter = env.solver().lra().getMaximalIterationCount();
}
// start the iterations
ValueType result = storm::utility::zero<ValueType>();
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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
void LraViHelper<ValueType, ComponentType, TransitionsType>::initializeNewValues(ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector<ValueType> 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<ValueType>() / _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<ValueType>());
_Tsx2 = _Tsx1;
if (_hasInstantStates) {
// Set-up vectors for storing intermediate results for instant states.
_Isx.resize(_IsTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
_Isb = _IsChoiceValues;
}
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
void LraViHelper<ValueType, ComponentType, TransitionsType>::prepareSolversAndMultipliers(const Environment& env, storm::solver::OptimizationDirection const* dir) {
_TsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _TsTransitions);
if (_hasInstantStates) {
if (_IsTransitions.getNonzeroEntryCount() > 0) {
// Set-up a solver for transitions within instant states
_IsSolverEnv = std::make_unique<storm::Environment>(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<ValueType> 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<ValueType> 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<ValueType> 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<ValueType>().create(env, _TsToIsTransitions);
_IsToTsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _IsToTsTransitions);
}
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
void LraViHelper<ValueType, ComponentType, TransitionsType>::setInputModelChoices(std::vector<uint64_t>& choices, std::vector<uint64_t> 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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
void LraViHelper<ValueType, ComponentType, TransitionsType>::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const* dir, std::vector<uint64_t>* 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<uint64_t> 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<uint64_t> 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 ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
typename LraViHelper<ValueType, ComponentType, TransitionsType>::ConvergenceCheckResult LraViHelper<ValueType, ComponentType, TransitionsType>::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<ValueType>() };
// Now check whether the currently produced results are precise enough
STORM_LOG_ASSERT(threshold > storm::utility::zero<ValueType>(), "Did not expect a non-positive threshold.");
auto x1It = xOld().begin();
auto x1Ite = xOld().end();
auto x2It = xNew().begin();
ValueType maxDiff = (*x2It - *x1It);
ValueType minDiff = maxDiff;
// The difference between maxDiff and minDiff is zero at this point. Thus, it doesn't make sense to check the threshold now.
for (++x1It, ++x2It; x1It != x1Ite; ++x1It, ++x2It) {
ValueType diff = (*x2It - *x1It);
// Potentially update maxDiff or minDiff
bool skipCheck = false;
if (maxDiff < diff) {
maxDiff = diff;
} else if (minDiff > diff) {
minDiff = diff;
} else {
skipCheck = true;
}
// Check convergence
if (!skipCheck && (maxDiff - minDiff) > (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<ValueType>(2.0));
// "Undo" the scaling of the values
res.currentValue = avgDiff * _uniformizationRate;
return res;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
void LraViHelper<ValueType, ComponentType, TransitionsType>::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<ValueType, ValueType>(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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
bool LraViHelper<ValueType, ComponentType, TransitionsType>::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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
std::vector<ValueType>& LraViHelper<ValueType, ComponentType, TransitionsType>::xNew() {
return _Tsx1IsCurrent ? _Tsx1 : _Tsx2;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
std::vector<ValueType> const& LraViHelper<ValueType, ComponentType, TransitionsType>::xNew() const {
return _Tsx1IsCurrent ? _Tsx1 : _Tsx2;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
std::vector<ValueType>& LraViHelper<ValueType, ComponentType, TransitionsType>::xOld() {
return _Tsx1IsCurrent ? _Tsx2 : _Tsx1;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
std::vector<ValueType> const& LraViHelper<ValueType, ComponentType, TransitionsType>::xOld() const {
return _Tsx1IsCurrent ? _Tsx2 : _Tsx1;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
bool LraViHelper<ValueType, ComponentType, TransitionsType>::nondetTs() const {
return TransitionsType == LraViTransitionsType::NondetTsNoIs;
}
template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
bool LraViHelper<ValueType, ComponentType, TransitionsType>::nondetIs() const {
return TransitionsType == LraViTransitionsType::DetTsNondetIs;
}
template class LraViHelper<double, storm::storage::MaximalEndComponent, LraViTransitionsType::NondetTsNoIs>;
template class LraViHelper<storm::RationalNumber, storm::storage::MaximalEndComponent, LraViTransitionsType::NondetTsNoIs>;
template class LraViHelper<double, storm::storage::MaximalEndComponent, LraViTransitionsType::DetTsNondetIs>;
template class LraViHelper<storm::RationalNumber, storm::storage::MaximalEndComponent, LraViTransitionsType::DetTsNondetIs>;
}
}
}
}

151
src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h
View File

@ -0,0 +1,151 @@
#pragma once
#include "storm/storage/SparseMatrix.h"
namespace storm {
class Environment;
namespace solver {
template<typename ValueType> class LinearEquationSolver;
template<typename ValueType> class MinMaxLinearEquationSolver;
template<typename ValueType> 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 <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
class LraViHelper {
public:
/// Function mapping from indices to values
typedef std::function<ValueType(uint64_t)> 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<ValueType> const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates = nullptr, std::vector<ValueType> 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<ValueType> const* exitRates = nullptr, storm::solver::OptimizationDirection const* dir = nullptr, std::vector<uint64_t>* 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<ValueType> 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<uint64_t>* 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<uint64_t>& choices, std::vector<uint64_t> 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<ValueType>& xNew();
std::vector<ValueType> const& xNew() const;
/// The result for timed states of the previous iteration
std::vector<ValueType>& xOld();
std::vector<ValueType> 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<ValueType> const& _transitionMatrix;
storm::storage::BitVector const* _timedStates; // e.g. Markovian states of a Markov automaton.
bool _hasInstantStates;
ValueType _uniformizationRate;
storm::storage::SparseMatrix<ValueType> _TsTransitions, _TsToIsTransitions, _IsTransitions, _IsToTsTransitions;
std::vector<ValueType> _Tsx1, _Tsx2, _TsChoiceValues;
bool _Tsx1IsCurrent;
std::vector<ValueType> _Isx, _Isb, _IsChoiceValues;
std::unique_ptr<storm::solver::Multiplier<ValueType>> _TsMultiplier, _TsToIsMultiplier, _IsToTsMultiplier;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> _NondetIsSolver;
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> _DetIsSolver;
std::unique_ptr<storm::Environment> _IsSolverEnv;
};
}
}
}
}
Write Preview
Loading…
Cancel
Save