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

NondeterministicLraHelper: Removed old ViHelper code.

tempestpy_adaptions
Tim Quatmann 4 years ago
parent
7e65e797fa
commit
f113ac7187
1 changed files with 18 additions and 533 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. 551
      src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp

551
src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp
View File

@ -1,10 +1,12 @@
#include "SparseNondeterministicInfiniteHorizonHelper.h"
#include "storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h"
#include "storm/solver/MinMaxLinearEquationSolver.h"
#include "storm/solver/LinearEquationSolver.h"
#include "storm/solver/Multiplier.h"
#include "storm/solver/LpSolver.h"
#include "storm/utility/graph.h"
#include "storm/utility/SignalHandler.h"
#include "storm/utility/solver.h"
#include "storm/utility/vector.h"
@ -84,9 +86,11 @@ namespace storm {
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.
underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
storm::RationalNumber newPrecision = env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2);
underlyingSolverEnvironment.solver().minMax().setPrecision(newPrecision);
underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion());
underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
underlyingSolverEnvironment.solver().setLinearEquationSolverPrecision(newPrecision, env.solver().lra().getRelativeTerminationCriterion());
underlyingSolverEnvironment.solver().lra().setPrecision(newPrecision);
}
// If requested, allocate memory for the choices made
@ -204,545 +208,26 @@ namespace storm {
}
}
/*!
* Abstract helper class that performs a single iteration of the value iteration method
*/
template <typename ValueType>
class LraViHelper {
public:
LraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix<ValueType> const& transitionMatrix) : _mec(mec), _transitionMatrix(transitionMatrix) {
// Intentionally left empty
}
virtual ~LraViHelper() = default;
/*!
* performs a single iteration step.
* If a choices vector is 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
* @return the current estimate of the LRA value
*/
virtual void iterate(Environment const& env, storm::solver::OptimizationDirection const& dir, std::vector<uint64_t>* choices = nullptr) = 0;
struct ConvergenceCheckResult {
bool isPrecisionAchieved;
ValueType currentValue;
};
/*!
* Checks whether the curently computed value achieves the desired precision
*/
virtual ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) = 0;
/*!
* Must be called between two calls of iterate.
*/
virtual void prepareNextIteration(Environment const& env) = 0;
protected:
/*!
*
* @param xPrevious the 'old' values
* @param xCurrent the 'new' values
* @param threshold the threshold
* @param relative whether the relative difference should be considered
* @return The first component is true if the (relative) difference between the maximal and the minimal entry-wise change of the two value vectors is below or equal to the provided threshold.
* In this case, the second component is the average of the maximal and the minimal change.
* If the threshold is exceeded, the computation is aborted early and the second component is only an approximation of the averages.
*/
std::pair<bool, ValueType> checkMinMaxDiffBelowThreshold(std::vector<ValueType> const& xPrevious, std::vector<ValueType> const& xCurrent, ValueType const& threshold, bool relative) const {
STORM_LOG_ASSERT(xPrevious.size() == xCurrent.size(), "Unexpected Dimension Mismatch");
STORM_LOG_ASSERT(threshold > storm::utility::zero<ValueType>(), "Did not expect a non-positive threshold.");
auto x1It = xPrevious.begin();
auto x1Ite = xPrevious.end();
auto x2It = xCurrent.begin();
ValueType maxDiff = (*x2It - *x1It);
ValueType minDiff = maxDiff;
bool result = true;
// 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)) {
result = false;
break;
}
}
ValueType avgDiff = (maxDiff + minDiff) / (storm::utility::convertNumber<ValueType>(2.0));
return {result, avgDiff};
}
storm::storage::MaximalEndComponent const& _mec;
storm::storage::SparseMatrix<ValueType> const& _transitionMatrix;
};
/*!
* Abstract helper class that performs a single iteration of the value iteration method for MDP
* @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
*/
template <typename ValueType>
class MdpLraViHelper : public LraViHelper<ValueType> {
public:
MdpLraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::function<ValueType(uint64_t stateIndex)> const& stateRewardsGetter, std::function<ValueType(uint64_t globalChoiceIndex)> const& actionRewardsGetter, ValueType const& aperiodicFactor) : LraViHelper<ValueType>(mec, transitionMatrix), _x1(mec.size(), storm::utility::zero<ValueType>()), _x2(_x1), _x1IsCurrent(true) {
// We add a selfloop to each state (which is necessary for convergence)
// Very roughly, this selfloop avoids that the values can flip around like this: [1, 0] -> [0, 1] -> [1, 0] -> ...
ValueType selfLoopProb = aperiodicFactor;
// Introducing the selfloop also requires the rewards to be scaled by the following factor.
_scalingFactor = storm::utility::one<ValueType>() - selfLoopProb;
uint64_t numMecStates = this->_mec.size();
boost::container::flat_map<uint64_t, uint64_t> toSubModelStateMapping;
toSubModelStateMapping.reserve(numMecStates);
uint64_t currState = 0;
uint64_t numMecChoices = 0;
for (auto const& stateChoices : this->_mec) {
toSubModelStateMapping.emplace(stateChoices.first, currState);
++currState;
numMecChoices += stateChoices.second.size();
}
assert(currState == numMecStates);
// Get a transition matrix that only considers the states and choices within the MEC
storm::storage::SparseMatrixBuilder<ValueType> mecTransitionBuilder(numMecChoices, numMecStates, 0, true, true, numMecStates);
_choiceValues.reserve(numMecChoices);
uint64_t currRow = 0;
for (auto const& stateChoices : this->_mec) {
auto const& mecState = stateChoices.first;
auto const& mecChoices = stateChoices.second;
mecTransitionBuilder.newRowGroup(currRow);
for (auto const& choice : mecChoices) {
bool insertedDiagElement = false;
for (auto const& entry : this->_transitionMatrix.getRow(choice)) {
uint64_t column = toSubModelStateMapping[entry.getColumn()];
if (!insertedDiagElement && entry.getColumn() > mecState) {
mecTransitionBuilder.addNextValue(currRow, toSubModelStateMapping[mecState], selfLoopProb);
insertedDiagElement = true;
}
if (!insertedDiagElement && entry.getColumn() == mecState) {
mecTransitionBuilder.addNextValue(currRow, column, selfLoopProb + _scalingFactor * entry.getValue());
insertedDiagElement = true;
} else {
mecTransitionBuilder.addNextValue(currRow, column, _scalingFactor * entry.getValue());
}
}
if (!insertedDiagElement) {
mecTransitionBuilder.addNextValue(currRow, toSubModelStateMapping[mecState], selfLoopProb);
}
// Compute the rewards obtained for this choice
_choiceValues.push_back(_scalingFactor * (stateRewardsGetter(mecState) + actionRewardsGetter(choice)));
++currRow;
}
}
_mecTransitions = mecTransitionBuilder.build();
STORM_LOG_ASSERT(_mecTransitions.isProbabilistic(), "The MEC-Matrix is not probabilistic.");
STORM_LOG_ASSERT(_mecTransitions.getRowGroupCount() == _x1.size(), "Unexpected size mismatch for created matrix.");
STORM_LOG_ASSERT(_x1.size() == _x2.size(), "Unexpected size mismatch for created matrix.");
}
virtual void iterate(Environment const& env, storm::solver::OptimizationDirection const& dir, std::vector<uint64_t>* choices = nullptr) override {
// Initialize a multipler if it does not exist, yet
if (!_multiplier) {
_multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _mecTransitions);
}
if (choices == nullptr) {
// Perform a simple matrix-vector multiplication
_multiplier->multiplyAndReduce(env, dir, xCurrent(), &_choiceValues, xPrevious());
} else {
// Perform a simple matrix-vector multiplication but also keep track of the choices within the _mecTransitions
std::vector<uint64_t> mecChoices(_mecTransitions.getRowGroupCount());
_multiplier->multiplyAndReduce(env, dir, xCurrent(), &_choiceValues, xPrevious(), &mecChoices);
// Transform the local choices (within this mec) to global indices
uint64_t mecState = 0;
for (auto const& stateChoices : this->_mec) {
uint64_t mecChoice = mecChoices[mecState];
STORM_LOG_ASSERT(mecChoice < stateChoices.second.size(), "The selected choice does not seem to exist.");
uint64_t globalChoiceIndex = *(stateChoices.second.begin() + mecChoice);
(*choices)[stateChoices.first] = globalChoiceIndex - this->_transitionMatrix.getRowGroupIndices()[stateChoices.first];
++mecState;
}
}
// Swap current and previous x vectors
_x1IsCurrent = !_x1IsCurrent;
}
virtual typename LraViHelper<ValueType>::ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) override {
typename LraViHelper<ValueType>::ConvergenceCheckResult res;
std::tie(res.isPrecisionAchieved, res.currentValue) = this->checkMinMaxDiffBelowThreshold(xPrevious(), xCurrent(), precision, relative);
res.currentValue /= _scalingFactor; // "Undo" the scaling of the rewards
return res;
}
virtual void prepareNextIteration(Environment const&) override {
// To avoid large (and numerically unstable) x-values, we substract a reference value.
ValueType referenceValue = xCurrent().front();
storm::utility::vector::applyPointwise<ValueType, ValueType>(xCurrent(), xCurrent(), [&referenceValue] (ValueType const& x_i) -> ValueType { return x_i - referenceValue; });
}
private:
std::vector<ValueType>& xCurrent() {
return _x1IsCurrent ? _x1 : _x2;
}
std::vector<ValueType>& xPrevious() {
return _x1IsCurrent ? _x2 : _x1;
}
storm::storage::SparseMatrix<ValueType> _mecTransitions;
std::vector<ValueType> _x1, _x2, _choiceValues;
bool _x1IsCurrent;
std::unique_ptr<storm::solver::Multiplier<ValueType>> _multiplier;
ValueType _scalingFactor;
};
/*!
* Abstract helper class that performs a single iteration of the value iteration method for MA
* @see Butkova, Wimmer, Hermanns: Long-Run Rewards for Markov Automata (TACAS'17), https://doi.org/10.1007/978-3-662-54580-5_11
*/
template <typename ValueType>
class MaLraViHelper : public LraViHelper<ValueType> {
public:
MaLraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector<ValueType> const& exitRates, std::function<ValueType(uint64_t stateIndex)> const& stateRewardsGetter, std::function<ValueType(uint64_t globalChoiceIndex)> const& actionRewardsGetter, ValueType const& aperiodicFactor) : LraViHelper<ValueType>(mec, transitionMatrix), _markovianStates(markovianStates), _Msx1IsCurrent(false) {
// Run through the Mec and collect some data:
// We consider two submodels, one consisting of the Markovian MEC states and one consisting of the probabilistic MEC states.
// 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 numPsSubModelStates(0), numPsSubModelChoices(0);
uint64_t numMsSubModelStates(0); // The number of choices coincide
// We will need to uniformize the Markovian 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 MEC.
_uniformizationRate = storm::utility::zero<ValueType>();
// Now run over the MEC and collect the required data.
for (auto const& stateChoices : this->_mec) {
uint64_t const& mecState = stateChoices.first;
if (_markovianStates.get(mecState)) {
toSubModelStateMapping.emplace(mecState, numMsSubModelStates);
++numMsSubModelStates;
STORM_LOG_ASSERT(stateChoices.second.size() == 1, "Markovian state has multiple MEC choices.");
_uniformizationRate = std::max(_uniformizationRate, exitRates[mecState]);
} else {
toSubModelStateMapping.emplace(mecState, numPsSubModelStates);
++numPsSubModelStates;
numPsSubModelChoices += stateChoices.second.size();
}
}
assert(numPsSubModelStates + numMsSubModelStates == mec.size());
STORM_LOG_THROW(numMsSubModelStates > 0, storm::exceptions::InvalidOperationException, "Markov Automaton has Zeno behavior. Computation of Long Run Average values not supported.");
_hasProbabilisticStates = numPsSubModelStates > 0;
// We make sure that every Markovian state gets a selfloop to make the model aperiodic
_uniformizationRate *= storm::utility::one<ValueType>() + aperiodicFactor;
// Now build the Markovian and the Probabilistic submodels.
// In addition, we also need the transitions between the two.
storm::storage::SparseMatrixBuilder<ValueType> msTransitionsBuilder(numMsSubModelStates, numMsSubModelStates);
_MsChoiceValues.reserve(numMsSubModelStates);
storm::storage::SparseMatrixBuilder<ValueType> msToPsTransitionsBuilder, psTransitionsBuilder, psToMsTransitionsBuilder;
if (_hasProbabilisticStates) {
msToPsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numMsSubModelStates, numPsSubModelStates);
psTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numPsSubModelChoices, numPsSubModelStates, 0, true, true, numPsSubModelStates);
psToMsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numPsSubModelChoices, numMsSubModelStates, 0, true, true, numPsSubModelStates);
_PsChoiceValues.reserve(numPsSubModelChoices);
}
uint64_t currMsRow = 0;
uint64_t currPsRow = 0;
for (auto const& stateChoices : this->_mec) {
uint64_t const& mecState = stateChoices.first;
auto const& mecChoices = stateChoices.second;
if (!_hasProbabilisticStates || _markovianStates.get(mecState)) {
// The currently processed state is Markovian.
// We need to uniformize!
ValueType uniformizationFactor = exitRates[mecState] / _uniformizationRate;
ValueType selfLoopProb = storm::utility::one<ValueType>() - uniformizationFactor;
STORM_LOG_ASSERT(mecChoices.size() == 1, "Unexpected number of choices at Markovian state.");
for (auto const& mecChoice : mecChoices) {
bool insertedDiagElement = false;
for (auto const& entry : this->_transitionMatrix.getRow(mecChoice)) {
uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()];
if (!_hasProbabilisticStates || _markovianStates.get(entry.getColumn())) {
// We have a transition from a Markovian state to a Markovian state
STORM_LOG_ASSERT(subModelColumn < numMsSubModelStates, "Invalid state for Markovian submodel");
if (!insertedDiagElement && subModelColumn > currMsRow) {
// We passed the diagonal entry, so add it now before moving on to the next entry
msTransitionsBuilder.addNextValue(currMsRow, currMsRow, selfLoopProb);
insertedDiagElement = true;
}
if (!insertedDiagElement && subModelColumn == currMsRow) {
// The current entry is the diagonal (selfloop) entry
msTransitionsBuilder.addNextValue(currMsRow, subModelColumn, selfLoopProb + uniformizationFactor * entry.getValue());
insertedDiagElement = true;
} else {
// The diagonal element either has been inserted already or still lies in front
msTransitionsBuilder.addNextValue(currMsRow, subModelColumn, uniformizationFactor * entry.getValue());
}
} else {
// We have a transition from a Markovian to a probabilistic state
STORM_LOG_ASSERT(subModelColumn < numPsSubModelStates, "Invalid state for probabilistic submodel");
msToPsTransitionsBuilder.addNextValue(currMsRow, subModelColumn, uniformizationFactor * entry.getValue());
}
}
// If the diagonal entry for the MS matrix still has not been set, we do that now
if (!insertedDiagElement) {
msTransitionsBuilder.addNextValue(currMsRow, currMsRow, selfLoopProb);
}
// Compute the rewards obtained for this choice.
_MsChoiceValues.push_back(stateRewardsGetter(mecState) / _uniformizationRate + actionRewardsGetter(mecChoice) * exitRates[mecState] / _uniformizationRate);
++currMsRow;
}
} else {
// The currently processed state is probabilistic
psTransitionsBuilder.newRowGroup(currPsRow);
psToMsTransitionsBuilder.newRowGroup(currPsRow);
for (auto const& mecChoice : mecChoices) {
for (auto const& entry : this->_transitionMatrix.getRow(mecChoice)) {
uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()];
if (_markovianStates.get(entry.getColumn())) {
// We have a transition from a probabilistic state to a Markovian state
STORM_LOG_ASSERT(subModelColumn < numMsSubModelStates, "Invalid state for Markovian submodel");
psToMsTransitionsBuilder.addNextValue(currPsRow, subModelColumn, entry.getValue());
} else {
// We have a transition from a probabilistic to a probabilistic state
STORM_LOG_ASSERT(subModelColumn < numPsSubModelStates, "Invalid state for probabilistic submodel");
psTransitionsBuilder.addNextValue(currPsRow, subModelColumn, entry.getValue());
}
}
// Compute the rewards obtained for this choice.
// State rewards do not count here since no time passes in probabilistic states.
_PsChoiceValues.push_back(actionRewardsGetter(mecChoice));
++currPsRow;
}
}
}
_MsTransitions = msTransitionsBuilder.build();
if (_hasProbabilisticStates) {
_MsToPsTransitions = msToPsTransitionsBuilder.build();
_PsTransitions = psTransitionsBuilder.build();
_PsToMsTransitions = psToMsTransitionsBuilder.build();
}
}
void initializeIterations(Environment const& env, storm::solver::OptimizationDirection const& dir) {
_Msx1.resize(_MsTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
_Msx2 = _Msx1;
_MsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _MsTransitions);
if (_hasProbabilisticStates) {
if (_PsTransitions.getNonzeroEntryCount() > 0) {
// Set-up a solver for transitions within PS states
_PsSolverEnv = env;
if (env.solver().isForceSoundness()) {
// To get correct results, the inner equation systems are solved exactly.
// TODO investigate how an error would propagate
_PsSolverEnv.solver().setForceExact(true);
}
storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> factory;
bool isAcyclic = !storm::utility::graph::hasCycle(_PsTransitions);
if (isAcyclic) {
STORM_LOG_INFO("Probabilistic transitions are acyclic.");
_PsSolverEnv.solver().minMax().setMethod(storm::solver::MinMaxMethod::Acyclic);
}
_PsSolver = factory.create(_PsSolverEnv, _PsTransitions);
_PsSolver->setHasUniqueSolution(true); // Assume non-zeno MA
_PsSolver->setHasNoEndComponents(true); // assume non-zeno MA
_PsSolver->setCachingEnabled(true);
_PsSolver->setRequirementsChecked(true);
auto req = _PsSolver->getRequirements(_PsSolverEnv, 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 checked.");
}
// Set up multipliers for transitions connecting Markovian and probabilistic states
_MsToPsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _MsToPsTransitions);
_PsToMsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _PsToMsTransitions);
// Set-up vectors for storing intermediate results for PS states.
_Psx.resize(_PsTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
_Psb = _PsChoiceValues;
}
}
void setInputModelChoices(std::vector<uint64_t>& choices, std::vector<uint64_t> const& localMecChoices, bool setChoiceZeroToMarkovianStates) {
// Transform the local choices (within this mec) to choice indices for the input model
uint64_t mecState = 0;
for (auto const& stateChoices : this->_mec) {
if (setChoiceZeroToMarkovianStates && _markovianStates.get(stateChoices.first)) {
choices[stateChoices.first] = 0;
} else {
uint64_t mecChoice = localMecChoices[mecState];
STORM_LOG_ASSERT(mecChoice < stateChoices.second.size(), "The selected choice does not seem to exist.");
uint64_t globalChoiceIndex = *(stateChoices.second.begin() + mecChoice);
choices[stateChoices.first] = globalChoiceIndex - this->_transitionMatrix.getRowGroupIndices()[stateChoices.first];
++mecState;
}
}
STORM_LOG_ASSERT(mecState == localMecChoices.size(), "Did not traverse all mec states.");
}
virtual void iterate(Environment const& env, storm::solver::OptimizationDirection const& dir, std::vector<uint64_t>* choices = nullptr) override {
// Initialize value vectors, multiplers, and solver if this has not been done, yet
if (!_MsMultiplier) {
initializeIterations(env, dir);
}
// Compute new x values for the Markovian states
// Flip what is current and what is previous
_Msx1IsCurrent = !_Msx1IsCurrent;
// At this point, xPrevious() points to what has been computed in the previous call of iterate (initially, this is the 0-vector).
// The result of this computation will be stored in xCurrent()
// Compute the values obtained by a single uniformization step between Markovian states only
_MsMultiplier->multiply(env, xPrevious(), &_MsChoiceValues, xCurrent());
if (_hasProbabilisticStates) {
// Add the values obtained by taking a single uniformization step that leads to a Probabilistic state followed by arbitrarily many probabilistic steps.
// First compute the total values when taking arbitrarily many probabilistic transitions (in no time)
if (_PsSolver) {
// We might need to track the optimal choices.
if (choices == nullptr) {
_PsSolver->solveEquations(_PsSolverEnv, dir, _Psx, _Psb);
} else {
_PsSolver->setTrackScheduler();
_PsSolver->solveEquations(_PsSolverEnv, dir, _Psx, _Psb);
setInputModelChoices(*choices, _PsSolver->getSchedulerChoices(), true);
}
} else {
STORM_LOG_ASSERT(_PsTransitions.getNonzeroEntryCount() == 0, "If no solver was initialized, an empty matrix would have been expected.");
if (choices == nullptr) {
storm::utility::vector::reduceVectorMinOrMax(dir, _Psb, _Psx, _PsTransitions.getRowGroupIndices());
} else {
std::vector<uint64_t> psMecChoices(_PsTransitions.getRowGroupCount());
storm::utility::vector::reduceVectorMinOrMax(dir, _Psb, _Psx, _PsTransitions.getRowGroupIndices(), &psMecChoices);
setInputModelChoices(*choices, _PsSolver->getSchedulerChoices(), true);
}
}
// Now add the (weighted) values of the probabilistic states to the values of the Markovian states.
_MsToPsMultiplier->multiply(env, _Psx, &xCurrent(), xCurrent());
}
}
virtual typename LraViHelper<ValueType>::ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) override {
typename LraViHelper<ValueType>::ConvergenceCheckResult res;
// 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);
std::tie(res.isPrecisionAchieved, res.currentValue) = this->checkMinMaxDiffBelowThreshold(xPrevious(), xCurrent(), threshold, relative);
res.currentValue *= _uniformizationRate; // "Undo" the scaling of the values
return res;
}
virtual void prepareNextIteration(Environment const& env) override {
// To avoid large (and numerically unstable) x-values, we substract a reference value.
ValueType referenceValue = xCurrent().front();
storm::utility::vector::applyPointwise<ValueType, ValueType>(xCurrent(), xCurrent(), [&referenceValue] (ValueType const& x_i) -> ValueType { return x_i - referenceValue; });
if (_hasProbabilisticStates) {
// Update the RHS of the equation system for the probabilistic states by taking the new values of Markovian states into account.
_PsToMsMultiplier->multiply(env, xCurrent(), &_PsChoiceValues, _Psb);
}
}
private:
std::vector<ValueType>& xCurrent() {
return _Msx1IsCurrent ? _Msx1 : _Msx2;
}
std::vector<ValueType>& xPrevious() {
return _Msx1IsCurrent ? _Msx2 : _Msx1;
}
storm::storage::BitVector const& _markovianStates;
bool _hasProbabilisticStates;
ValueType _uniformizationRate;
storm::storage::SparseMatrix<ValueType> _MsTransitions, _MsToPsTransitions, _PsTransitions, _PsToMsTransitions;
std::vector<ValueType> _Msx1, _Msx2, _MsChoiceValues;
bool _Msx1IsCurrent;
std::vector<ValueType> _Psx, _Psb, _PsChoiceValues;
std::unique_ptr<storm::solver::Multiplier<ValueType>> _MsMultiplier, _MsToPsMultiplier, _PsToMsMultiplier;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> _PsSolver;
Environment _PsSolverEnv;
};
template <typename ValueType>
ValueType SparseNondeterministicInfiniteHorizonHelper<ValueType>::computeLraForMecVi(Environment const& env, std::function<ValueType(uint64_t stateIndex)> const& stateRewardsGetter, std::function<ValueType(uint64_t globalChoiceIndex)> const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) {
// Collect some parameters of the computation.
// Collect some parameters of the computation
ValueType aperiodicFactor = storm::utility::convertNumber<ValueType>(env.solver().lra().getAperiodicFactor());
ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision()) / aperiodicFactor;
bool relative = env.solver().lra().getRelativeTerminationCriterion();
boost::optional<uint64_t> maxIter;
if (env.solver().lra().isMaximalIterationCountSet()) {
maxIter = env.solver().lra().getMaximalIterationCount();
std::vector<uint64_t>* optimalChoices = nullptr;
if (isProduceSchedulerSet()) {
optimalChoices = &_producedOptimalChoices.get();
}
auto dir = this->getOptimizationDirection();
// Create an object for the iterations
std::shared_ptr<LraViHelper<ValueType>> iterationHelper;
// Now create a helper and perform the algorithm
if (isContinuousTime()) {
iterationHelper = std::make_shared<MaLraViHelper<ValueType>>(mec, _transitionMatrix, *_markovianStates, *_exitRates, stateRewardsGetter, actionRewardsGetter, aperiodicFactor);
} else {
iterationHelper = std::make_shared<MdpLraViHelper<ValueType>>(mec, _transitionMatrix, stateRewardsGetter, actionRewardsGetter, aperiodicFactor);
}
// start the iterations
ValueType result = storm::utility::zero<ValueType>();
uint64_t iter = 0;
while (!maxIter.is_initialized() || iter < maxIter.get()) {
++iter;
iterationHelper->iterate(env, dir);
// Check if we are done
auto convergenceCheckResult = iterationHelper->checkConvergence(relative, precision);
result = convergenceCheckResult.currentValue;
if (convergenceCheckResult.isPrecisionAchieved) {
break;
}
if (storm::utility::resources::isTerminate()) {
break;
}
iterationHelper->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.");
// We assume a Markov Automaton (with deterministic timed states and nondeterministic instant states)
storm::modelchecker::helper::internal::LraViHelper<ValueType, storm::storage::MaximalEndComponent, storm::modelchecker::helper::internal::LraViTransitionsType::DetTsNondetIs> viHelper(mec, _transitionMatrix, aperiodicFactor, _markovianStates, _exitRates);
return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, _exitRates, &this->getOptimizationDirection(), optimalChoices);
} else {
STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations.");
// We assume an MDP (with nondeterministic timed states and no instant states)
storm::modelchecker::helper::internal::LraViHelper<ValueType, storm::storage::MaximalEndComponent, storm::modelchecker::helper::internal::LraViTransitionsType::NondetTsNoIs> viHelper(mec, _transitionMatrix, aperiodicFactor);
return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, nullptr, &this->getOptimizationDirection(), optimalChoices);
}
if (isProduceSchedulerSet()) {
// We will be doing one more iteration step and track scheduler choices this time.
iterationHelper->prepareNextIteration(env);
iterationHelper->iterate(env, dir, &_producedOptimalChoices.get());
}
return result;
}
template <typename ValueType>

Write Preview
Loading…
Cancel
Save