diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp index 979337f5f..7a9cb003e 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.cpp @@ -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(2)); + storm::RationalNumber newPrecision = env.solver().lra().getPrecision() / storm::utility::convertNumber(2); + underlyingSolverEnvironment.solver().minMax().setPrecision(newPrecision); underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion()); - underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber(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 - class LraViHelper { - public: - LraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix 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* 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 checkMinMaxDiffBelowThreshold(std::vector const& xPrevious, std::vector 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(), "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(2.0)); - return {result, avgDiff}; - } - - storm::storage::MaximalEndComponent const& _mec; - storm::storage::SparseMatrix 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 - class MdpLraViHelper : public LraViHelper { - public: - - MdpLraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix const& transitionMatrix, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, ValueType const& aperiodicFactor) : LraViHelper(mec, transitionMatrix), _x1(mec.size(), storm::utility::zero()), _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() - selfLoopProb; - - uint64_t numMecStates = this->_mec.size(); - boost::container::flat_map 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 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* choices = nullptr) override { - // Initialize a multipler if it does not exist, yet - if (!_multiplier) { - _multiplier = storm::solver::MultiplierFactory().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 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::ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) override { - typename LraViHelper::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(xCurrent(), xCurrent(), [&referenceValue] (ValueType const& x_i) -> ValueType { return x_i - referenceValue; }); - } - - private: - - std::vector& xCurrent() { - return _x1IsCurrent ? _x1 : _x2; - } - - std::vector& xPrevious() { - return _x1IsCurrent ? _x2 : _x1; - } - - storm::storage::SparseMatrix _mecTransitions; - std::vector _x1, _x2, _choiceValues; - bool _x1IsCurrent; - std::unique_ptr> _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 - class MaLraViHelper : public LraViHelper { - public: - - MaLraViHelper(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, std::vector const& exitRates, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, ValueType const& aperiodicFactor) : LraViHelper(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 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(); - // 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() + aperiodicFactor; - - // Now build the Markovian and the Probabilistic submodels. - // In addition, we also need the transitions between the two. - storm::storage::SparseMatrixBuilder msTransitionsBuilder(numMsSubModelStates, numMsSubModelStates); - _MsChoiceValues.reserve(numMsSubModelStates); - storm::storage::SparseMatrixBuilder msToPsTransitionsBuilder, psTransitionsBuilder, psToMsTransitionsBuilder; - if (_hasProbabilisticStates) { - msToPsTransitionsBuilder = storm::storage::SparseMatrixBuilder(numMsSubModelStates, numPsSubModelStates); - psTransitionsBuilder = storm::storage::SparseMatrixBuilder(numPsSubModelChoices, numPsSubModelStates, 0, true, true, numPsSubModelStates); - psToMsTransitionsBuilder = storm::storage::SparseMatrixBuilder(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() - 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()); - _Msx2 = _Msx1; - _MsMultiplier = storm::solver::MultiplierFactory().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 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().create(env, _MsToPsTransitions); - _PsToMsMultiplier = storm::solver::MultiplierFactory().create(env, _PsToMsTransitions); - - // Set-up vectors for storing intermediate results for PS states. - _Psx.resize(_PsTransitions.getRowGroupCount(), storm::utility::zero()); - _Psb = _PsChoiceValues; - } - - } - - void setInputModelChoices(std::vector& choices, std::vector 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* 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 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::ConvergenceCheckResult checkConvergence(bool relative, ValueType precision) override { - typename LraViHelper::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(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& xCurrent() { - return _Msx1IsCurrent ? _Msx1 : _Msx2; - } - - std::vector& xPrevious() { - return _Msx1IsCurrent ? _Msx2 : _Msx1; - } - - storm::storage::BitVector const& _markovianStates; - bool _hasProbabilisticStates; - ValueType _uniformizationRate; - storm::storage::SparseMatrix _MsTransitions, _MsToPsTransitions, _PsTransitions, _PsToMsTransitions; - std::vector _Msx1, _Msx2, _MsChoiceValues; - bool _Msx1IsCurrent; - std::vector _Psx, _Psb, _PsChoiceValues; - std::unique_ptr> _MsMultiplier, _MsToPsMultiplier, _PsToMsMultiplier; - std::unique_ptr> _PsSolver; - Environment _PsSolverEnv; - }; - template ValueType SparseNondeterministicInfiniteHorizonHelper::computeLraForMecVi(Environment const& env, std::function const& stateRewardsGetter, std::function const& actionRewardsGetter, storm::storage::MaximalEndComponent const& mec) { - // Collect some parameters of the computation. + // Collect some parameters of the computation ValueType aperiodicFactor = storm::utility::convertNumber(env.solver().lra().getAperiodicFactor()); - ValueType precision = storm::utility::convertNumber(env.solver().lra().getPrecision()) / aperiodicFactor; - bool relative = env.solver().lra().getRelativeTerminationCriterion(); - boost::optional maxIter; - if (env.solver().lra().isMaximalIterationCountSet()) { - maxIter = env.solver().lra().getMaximalIterationCount(); + std::vector* optimalChoices = nullptr; + if (isProduceSchedulerSet()) { + optimalChoices = &_producedOptimalChoices.get(); } - auto dir = this->getOptimizationDirection(); - // Create an object for the iterations - std::shared_ptr> iterationHelper; + // Now create a helper and perform the algorithm if (isContinuousTime()) { - iterationHelper = std::make_shared>(mec, _transitionMatrix, *_markovianStates, *_exitRates, stateRewardsGetter, actionRewardsGetter, aperiodicFactor); - } else { - iterationHelper = std::make_shared>(mec, _transitionMatrix, stateRewardsGetter, actionRewardsGetter, aperiodicFactor); - } - - // start the iterations - ValueType result = storm::utility::zero(); - 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 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 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