From bcd193dd5711445a00648d84931ff0e5669eb663 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Mon, 2 Dec 2019 13:53:23 +0100 Subject: [PATCH] Implemented Value iteration based LRA computation for CTMCs. --- .../csl/SparseCtmcCslModelChecker.cpp | 6 +- .../csl/helper/HybridCtmcCslHelper.cpp | 9 +- .../csl/helper/SparseCtmcCslHelper.cpp | 493 +++++++++++------- .../csl/helper/SparseCtmcCslHelper.h | 19 +- src/storm/storage/SparseMatrix.cpp | 46 +- src/storm/storage/SparseMatrix.h | 13 + 6 files changed, 366 insertions(+), 220 deletions(-) diff --git a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 89a2c9e03..05555b125 100644 --- a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -132,16 +132,14 @@ namespace storm { std::unique_ptr subResultPointer = this->check(env, stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - storm::storage::SparseMatrix probabilityMatrix = storm::modelchecker::helper::SparseCtmcCslHelper::computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); - std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal(this->getModel(), checkTask), probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector()); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), &this->getModel().getExitRateVector()); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template std::unique_ptr SparseCtmcCslModelChecker::computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { - storm::storage::SparseMatrix probabilityMatrix = storm::modelchecker::helper::SparseCtmcCslHelper::computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask); - std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal(this->getModel(), checkTask), probabilityMatrix, rewardModel.get(), &this->getModel().getExitRateVector()); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), rewardModel.get(), &this->getModel().getExitRateVector()); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } diff --git a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp index d76a28884..3555c7216 100644 --- a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp @@ -339,19 +339,18 @@ namespace storm { template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates) { - storm::dd::Add probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector); storm::utility::Stopwatch conversionWatch(true); // Create ODD for the translation. storm::dd::Odd odd = model.getReachableStates().createOdd(); - storm::storage::SparseMatrix explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd); + storm::storage::SparseMatrix explicitRateMatrix = rateMatrix.toMatrix(odd, odd); std::vector explicitExitRateVector = exitRateVector.toVector(odd); conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal(), explicitProbabilityMatrix, psiStates.toVector(odd), &explicitExitRateVector); + std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal(), explicitRateMatrix, psiStates.toVector(odd), &explicitExitRateVector); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); } @@ -367,12 +366,12 @@ namespace storm { // Create ODD for the translation. storm::dd::Odd odd = model.getReachableStates().createOdd(); - storm::storage::SparseMatrix explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd); + storm::storage::SparseMatrix explicitRateMatrix = rateMatrix.toMatrix(odd, odd); std::vector explicitExitRateVector = exitRateVector.toVector(odd); conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal(), explicitProbabilityMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector); + std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal(), explicitRateMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); } diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index ed87a9e1a..2125ee677 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -25,6 +25,7 @@ #include "storm/exceptions/InvalidPropertyException.h" #include "storm/exceptions/FormatUnsupportedBySolverException.h" #include "storm/exceptions/UncheckedRequirementException.h" +#include "storm/exceptions/NotSupportedException.h" namespace storm { namespace modelchecker { @@ -370,10 +371,10 @@ namespace storm { } template - std::vector SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector) { + std::vector SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector) { // If there are no goal states, we avoid the computation and directly return zero. - uint_fast64_t numberOfStates = probabilityMatrix.getRowCount(); + uint_fast64_t numberOfStates = rateMatrix.getRowCount(); if (psiStates.empty()) { return std::vector(numberOfStates, storm::utility::zero()); } @@ -386,7 +387,7 @@ namespace storm { ValueType zero = storm::utility::zero(); ValueType one = storm::utility::one(); - return computeLongRunAverages(env, std::move(goal), probabilityMatrix, + return computeLongRunAverages(env, std::move(goal), rateMatrix, [&zero, &one, &psiStates] (storm::storage::sparse::state_type const& state) -> ValueType { if (psiStates.get(state)) { return one; @@ -397,16 +398,33 @@ namespace storm { } template - std::vector SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, RewardModelType const& rewardModel, std::vector const* exitRateVector) { + std::vector SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, RewardModelType const& rewardModel, std::vector const* exitRateVector) { // Only compute the result if the model has a state-based reward model. STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - return computeLongRunAverageRewards(env, std::move(goal), probabilityMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, *exitRateVector), exitRateVector); + return computeLongRunAverages(env, std::move(goal), rateMatrix, + [&] (storm::storage::sparse::state_type const& state) -> ValueType { + ValueType result = rewardModel.hasStateRewards() ? rewardModel.getStateReward(state) : storm::utility::zero(); + if (rewardModel.hasStateActionRewards()) { + // State action rewards are multiplied with the exit rate r(s). Then, multiplying the reward with the expected time we stay at s (i.e. 1/r(s)) yields the original state reward + if (exitRateVector) { + result += rewardModel.getStateActionReward(state) * (*exitRateVector)[state]; + } else { + result += rewardModel.getStateActionReward(state); + } + } + if (rewardModel.hasTransitionRewards()) { + // Transition rewards are already multiplied with the rates + result += rateMatrix.getPointwiseProductRowSum(rewardModel.getTransitionRewardMatrix(), state); + } + return result; + }, + exitRateVector); } template - std::vector SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, std::vector const& stateRewardVector, std::vector const* exitRateVector) { - return computeLongRunAverages(env, std::move(goal), probabilityMatrix, + std::vector SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, std::vector const& stateRewardVector, std::vector const* exitRateVector) { + return computeLongRunAverages(env, std::move(goal), rateMatrix, [&stateRewardVector] (storm::storage::sparse::state_type const& state) -> ValueType { return stateRewardVector[state]; }, @@ -414,183 +432,49 @@ namespace storm { } template - std::vector SparseCtmcCslHelper::computeLongRunAverages(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, std::function const& valueGetter, std::vector const* exitRateVector){ - uint_fast64_t numberOfStates = probabilityMatrix.getRowCount(); - + std::vector SparseCtmcCslHelper::computeLongRunAverages(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, std::function const& valueGetter, std::vector const* exitRateVector){ + storm::storage::SparseMatrix probabilityMatrix; + if (exitRateVector) { + probabilityMatrix = computeProbabilityMatrix(rateMatrix, *exitRateVector); + } else { + probabilityMatrix = rateMatrix; + } + uint_fast64_t numberOfStates = rateMatrix.getRowCount(); + // Start by decomposing the CTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(probabilityMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs()); + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(rateMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs()); STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs."); - - // Get some data members for convenience. - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); - + // Prepare the vector holding the LRA values for each of the BSCCs. - std::vector bsccLra(bsccDecomposition.size(), zero); + std::vector bsccLra; + bsccLra.reserve(bsccDecomposition.size()); - // First we check which states are in BSCCs. + // Keep track of the maximal and minimal value occuring in one of the BSCCs + ValueType maxValue, minValue; storm::storage::BitVector statesInBsccs(numberOfStates); - storm::storage::BitVector firstStatesInBsccs(numberOfStates); - - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // Gather information for later use. - bool first = true; + auto backwardTransitions = rateMatrix.transpose(); + for (auto const& bscc : bsccDecomposition) { for (auto const& state : bscc) { statesInBsccs.set(state); - if (first) { - firstStatesInBsccs.set(state); - } - first = false; + } + bsccLra.push_back(computeLongRunAveragesForBscc(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector)); + if (bsccLra.size() == 1) { + maxValue = bsccLra.back(); + minValue = bsccLra.back(); + } else { + maxValue = std::max(bsccLra.back(), maxValue); + minValue = std::min(bsccLra.back(), minValue); } } - storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; + storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs."); - // Prepare a vector holding the index within all states that are in BSCCs for every state. - std::vector indexInStatesInBsccs; - - // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC. - std::vector stateToBsccIndexMap; - - if (!statesInBsccs.empty()) { - firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; - - // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. - storm::storage::SparseMatrix bsccEquationSystem = probabilityMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); - - // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a - // multiplication from the right by transposing the system. - bsccEquationSystem = bsccEquationSystem.transpose(false, true); - - // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states. - uint_fast64_t lastIndex = 0; - uint_fast64_t currentNumberOfSetBits = 0; - indexInStatesInBsccs.reserve(probabilityMatrix.getRowCount()); - for (auto index : statesInBsccs) { - while (lastIndex <= index) { - indexInStatesInBsccs.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } - - stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits()); - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - for (auto const& state : bscc) { - stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex; - } - } - - // Build a different system depending on the problem format of the equation solver. - // Check solver requirements. - storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; - auto requirements = linearEquationSolverFactory.getRequirements(env); - requirements.clearLowerBounds(); - requirements.clearUpperBounds(); - STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); - - bool fixedPointSystem = false; - if (linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem) { - fixedPointSystem = true; - } - - // Now build the final equation system matrix, the initial guess and the right-hand side in one go. - std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); - storm::storage::SparseMatrixBuilder builder; - for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { - - // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the - // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the - // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal. - if (firstStatesInBsccs.get(row)) { - uint_fast64_t requiredBscc = stateToBsccIndexMap[row]; - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc]; - - if (fixedPointSystem) { - for (auto const& state : bscc) { - if (row == indexInStatesInBsccs[state]) { - builder.addNextValue(row, indexInStatesInBsccs[state], zero); - } else { - builder.addNextValue(row, indexInStatesInBsccs[state], -one); - } - } - } else { - for (auto const& state : bscc) { - builder.addNextValue(row, indexInStatesInBsccs[state], one); - } - } - - bsccEquationSystemRightSide[row] = one; - } else { - // Otherwise, we copy the row, and subtract 1 from the diagonal (only for the equation solver format). - for (auto& entry : bsccEquationSystem.getRow(row)) { - if (fixedPointSystem || entry.getColumn() != row) { - builder.addNextValue(row, entry.getColumn(), entry.getValue()); - } else { - builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); - } - } - } - } - - // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC. - std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - for (auto const& state : bscc) { - bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); - } - } - - bsccEquationSystem = builder.build(); - { - std::unique_ptr> solver = linearEquationSolverFactory.create(env, std::move(bsccEquationSystem)); - solver->setLowerBound(storm::utility::zero()); - solver->setUpperBound(storm::utility::one()); - solver->solveEquations(env, bsccEquationSystemSolution, bsccEquationSystemRightSide); - } - - // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. - if (exitRateVector != nullptr) { - std::vector bsccTotalValue(bsccDecomposition.size(), zero); - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { - bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]); - } - - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { - bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]]; - } - } - - // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - for (auto const& state : bscc) { - bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += valueGetter(state) * bsccEquationSystemSolution[indexInStatesInBsccs[state]]; - } - } - - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); - } - } else { - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values - // directly (based on whether or not the contained state is a psi state). - bsccLra[bsccIndex] = valueGetter(*bscc.begin()); - } - - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); + std::vector stateToBsccMap(statesInBsccs.size(), -1); + for (uint64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + for (auto const& state : bsccDecomposition[bsccIndex]) { + stateToBsccMap[state] = bsccIndex; } } @@ -605,30 +489,44 @@ namespace storm { rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); for (auto state : statesNotInBsccs) { - ValueType reward = zero; - for (auto entry : probabilityMatrix.getRow(state)) { + ValueType reward = storm::utility::zero(); + for (auto entry : rateMatrix.getRow(state)) { if (statesInBsccs.get(entry.getColumn())) { - reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]]; + if (exitRateVector) { + reward += (entry.getValue() / (*exitRateVector)[state]) * bsccLra[stateToBsccMap[entry.getColumn()]]; + } else { + reward += entry.getValue() * bsccLra[stateToBsccMap[entry.getColumn()]]; + } } } rewardRightSide.push_back(reward); } - storm::storage::SparseMatrix rewardEquationSystemMatrix = probabilityMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); - rewardEquationSystemMatrix.convertToEquationSystem(); - - rewardSolution = std::vector(rewardEquationSystemMatrix.getColumnCount(), one); - - { - storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; - // Check solver requirements - auto requirements = linearEquationSolverFactory.getRequirements(env); - STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); - // Check whether we have the right input format for the solver. - STORM_LOG_THROW(linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem, storm::exceptions::FormatUnsupportedBySolverException, "The selected solver does not support the required format."); - std::unique_ptr> solver = linearEquationSolverFactory.create(env, std::move(rewardEquationSystemMatrix)); - solver->solveEquations(env, rewardSolution, rewardRightSide); + // Compute reachability rewards + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + bool isEqSysFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + storm::storage::SparseMatrix rewardEquationSystemMatrix = rateMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, isEqSysFormat); + if (exitRateVector) { + uint64_t localRow = 0; + for (auto const& globalRow : statesNotInBsccs) { + for (auto& entry : rewardEquationSystemMatrix.getRow(localRow)) { + entry.setValue(entry.getValue() / (*exitRateVector)[globalRow]); + } + ++localRow; + } + } + if (isEqSysFormat) { + rewardEquationSystemMatrix.convertToEquationSystem(); } + rewardSolution = std::vector(rewardEquationSystemMatrix.getColumnCount(), (maxValue + minValue) / storm::utility::convertNumber(2)); + // std::cout << rewardEquationSystemMatrix << std::endl; + // std::cout << storm::utility::vector::toString(rewardRightSide) << std::endl; + std::unique_ptr> solver = linearEquationSolverFactory.create(env, std::move(rewardEquationSystemMatrix)); + solver->setBounds(minValue, maxValue); + // Check solver requirements + auto requirements = solver->getRequirements(env); + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + solver->solveEquations(env, rewardSolution, rewardRightSide); } // Fill the result vector. @@ -653,6 +551,225 @@ namespace storm { return result; } + template + ValueType SparseCtmcCslHelper::computeLongRunAveragesForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector) { + std::cout << std::endl << "========== BSCC =======" << std::endl; + storm::storage::BitVector bsccStates(rateMatrix.getRowCount()); + for (auto const& s : bscc) { bsccStates.set(s); std::cout << s << ": " << valueGetter(s) << "\t";} + auto subm = rateMatrix.getSubmatrix(false,bsccStates,bsccStates); + std::cout << std::endl << subm << std::endl; + + std::cout << std::endl << "========== EQSYS =======" << std::endl; + + // Catch the trivial case for bsccs of size 1 + if (bscc.size() == 1) { + std::cout << std::endl << "========== Result = " << valueGetter(*bscc.begin()) << " =======" << std::endl; + + return valueGetter(*bscc.begin()); + } + + return computeLongRunAveragesForBsccVi(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector); + return computeLongRunAveragesForBsccEqSys(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector); + } + + template <> + storm::RationalFunction SparseCtmcCslHelper::computeLongRunAveragesForBsccVi(Environment const&, storm::storage::StronglyConnectedComponent const&, storm::storage::SparseMatrix const&, storm::storage::SparseMatrix const&, std::function const&, std::vector const*) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The requested Method for LRA computation is not supported for parametric models."); + } + + template + ValueType SparseCtmcCslHelper::computeLongRunAveragesForBsccVi(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector) { + + // Initialize data about the bscc + storm::storage::BitVector bsccStates(rateMatrix.getRowGroupCount(), false); + for (auto const& state : bscc) { + bsccStates.set(state); + } + + // Get the uniformization rate + ValueType uniformizationRate = storm::utility::one(); + if (exitRateVector) { + uniformizationRate = storm::utility::vector::max_if(*exitRateVector, bsccStates); + } + // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop. + // Hence, we increase the uniformization rate a little. + uniformizationRate += storm::utility::one(); // Todo: try other values such as *=1.01 + + // Get the transitions of the submodel + typename storm::storage::SparseMatrix bsccMatrix = rateMatrix.getSubmatrix(true, bsccStates, bsccStates, true); + + // Uniformize the transitions + uint64_t subState = 0; + for (auto state : bsccStates) { + for (auto& entry : bsccMatrix.getRow(subState)) { + if (entry.getColumn() == subState) { + if (exitRateVector) { + entry.setValue(storm::utility::one() + (entry.getValue() - (*exitRateVector)[state]) / uniformizationRate); + } else { + entry.setValue(storm::utility::one() + (entry.getValue() - storm::utility::one()) / uniformizationRate); + } + } else { + entry.setValue(entry.getValue() / uniformizationRate); + } + } + ++subState; + } + + // Compute the rewards obtained in a single uniformization step + std::vector markovianRewards; + markovianRewards.reserve(bsccMatrix.getRowCount()); + for (auto const& state : bsccStates) { + ValueType stateRewardScalingFactor = storm::utility::one() / uniformizationRate; + // ValueType actionRewardScalingFactor = exitRateVector[state] / uniformizationRate; // action rewards are currently not supported + markovianRewards.push_back(valueGetter(state) * stateRewardScalingFactor); + } + + // start the iterations + ValueType precision = storm::utility::convertNumber(1e-6) / uniformizationRate; // TODO env.solver().minMax().getPrecision()) / uniformizationRate; + bool relative = true; // TODO env.solver().minMax().getRelativeTerminationCriterion(); + std::vector v(bsccMatrix.getRowCount(), storm::utility::zero()); + std::vector w = v; + std::unique_ptr> solver; + while (true) { + // Compute the values for the markovian states. We also keep track of the maximal and minimal difference between two values (for convergence checking) + auto vIt = v.begin(); + uint64_t row = 0; + ValueType newValue = markovianRewards[row] + bsccMatrix.multiplyRowWithVector(row, w); + ValueType maxDiff = newValue - *vIt; + ValueType minDiff = maxDiff; + *vIt = newValue; + for (++vIt, ++row; row < bsccMatrix.getRowCount(); ++vIt, ++row) { + newValue = markovianRewards[row] + bsccMatrix.multiplyRowWithVector(row, w); + ValueType diff = newValue - *vIt; + maxDiff = std::max(maxDiff, diff); + minDiff = std::min(minDiff, diff); + *vIt = newValue; + } + + // Check for convergence + if ((maxDiff - minDiff) <= (relative ? (precision * (v.front() + minDiff)) : precision)) { + break; + } + + // update the rhs of the equation system + ValueType referenceValue = v.front(); + storm::utility::vector::applyPointwise(v, w, [&referenceValue] (ValueType const& v_i) -> ValueType { return v_i - referenceValue; }); + } + std::cout << std::endl << "========== VI Result = " << v.front() * uniformizationRate << " =======" << std::endl; + return v.front() * uniformizationRate; + } + + template + ValueType SparseCtmcCslHelper::computeLongRunAveragesForBsccEqSys(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector) { + + // Get a mapping from global state indices to local ones. + std::unordered_map toLocalIndexMap; + uint64_t localIndex = 0; + for (auto const& globalIndex : bscc) { + toLocalIndexMap[globalIndex] = localIndex; + ++localIndex; + } + + // Build the equation system matrix + + + // Build the equation system matrix for A[s,s] = R(s,s) - r(s) & A[s,s'] = R(s,s') (where s != s') + // x_0+...+x_n=1 & x*A=0 <=> x_0+...+x_n=1 & A^t*x=0 [ <=> 1-x_1+...+x_n=x_0 & (1-A^t)*x = x ] + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + storm::storage::SparseMatrixBuilder builder(bscc.size(), bscc.size()); + bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + // The first row asserts that the values sum up to one + uint64_t row = 0; + if (isEquationSystemFormat) { + for (uint64_t state = 0; state < bscc.size(); ++state) { + builder.addNextValue(row, state, storm::utility::one()); + } + } else { + for (uint64_t state = 1; state < bscc.size(); ++state) { + builder.addNextValue(row, state, -storm::utility::one()); + } + } + ++row; + // Build the equation system matrix + // We can skip the first state to make the equation system matrix square + auto stateIt = bscc.begin(); + for (++stateIt; stateIt != bscc.end(); ++stateIt) { + ValueType diagonalValue = (exitRateVector == nullptr) ? -storm::utility::one() : -(*exitRateVector)[*stateIt]; + if (!isEquationSystemFormat) { + diagonalValue = storm::utility::one() - diagonalValue; + } + bool insertedDiagonal = storm::utility::isZero(diagonalValue); + for (auto const& backwardsEntry : backwardTransitions.getRow(*stateIt)) { + auto localIndexIt = toLocalIndexMap.find(backwardsEntry.getColumn()); + if (localIndexIt != toLocalIndexMap.end()) { + ValueType val = backwardsEntry.getValue(); + if (!isEquationSystemFormat) { + val = -val; + } + uint64_t localIndex = localIndexIt->second; + if (!insertedDiagonal && localIndex == row) { + builder.addNextValue(row, localIndex, val + diagonalValue); + insertedDiagonal = true; + } else { + if (!insertedDiagonal && localIndex > row) { + builder.addNextValue(row, row, diagonalValue); + insertedDiagonal = true; + } + builder.addNextValue(row, localIndex, val); + } + } + } + if (!insertedDiagonal) { + builder.addNextValue(row, row, diagonalValue); + } + ++row; + } + + // Create a linear equation solver + auto matrix = builder.build(); + std::cout << matrix << std::endl; + auto solver = linearEquationSolverFactory.create(env, std::move(matrix)); + solver->setBounds(storm::utility::zero(), storm::utility::one()); + // Check solver requirements. + auto requirements = solver->getRequirements(env); + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + + std::vector bsccEquationSystemRightSide(bscc.size(), storm::utility::zero()); + bsccEquationSystemRightSide.front() = storm::utility::one(); + std::vector bsccEquationSystemSolution(bscc.size(), storm::utility::one() / storm::utility::convertNumber(bscc.size())); + std::cout << storm::utility::vector::toString(bsccEquationSystemRightSide) << std::endl; + solver->solveEquations(env, bsccEquationSystemSolution, bsccEquationSystemRightSide); + std::cout << storm::utility::vector::toString(bsccEquationSystemSolution) << std::endl; + // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. + if (false && exitRateVector != nullptr) { + ValueType totalValue = storm::utility::zero(); + auto solIt = bsccEquationSystemSolution.begin(); + for (auto const& globalState : bscc) { + totalValue += (*solIt) * (storm::utility::one() / (*exitRateVector)[globalState]); + ++solIt; + } + assert(solIt == bsccEquationSystemSolution.end()); + solIt = bsccEquationSystemSolution.begin(); + for (auto const& globalState : bscc) { + *solIt = ((*solIt) * (storm::utility::one() / (*exitRateVector)[globalState])) / totalValue; + ++solIt; + } + assert(solIt == bsccEquationSystemSolution.end()); + } + + // Calculate final LRA Value + ValueType result = storm::utility::zero(); + auto solIt = bsccEquationSystemSolution.begin(); + for (auto const& globalState : bscc) { + result += valueGetter(globalState) * (*solIt); + ++solIt; + } + assert(solIt == bsccEquationSystemSolution.end()); + std::cout << std::endl << "========== Result = " << result << " =======" << std::endl; + + return result; + } + template ::SupportsExponential, int>::type> std::vector SparseCtmcCslHelper::computeAllTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, double timeBound) { diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h index 4fb01d5ff..8a45ea75f 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h @@ -14,6 +14,10 @@ namespace storm { class Environment; + namespace storage { + class StronglyConnectedComponent; + } + namespace modelchecker { namespace helper { class SparseCtmcCslHelper { @@ -52,13 +56,13 @@ namespace storm { static std::vector computeTotalRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, RewardModelType const& rewardModel, bool qualitative); template - static std::vector computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector); + static std::vector computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector); template - static std::vector computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, RewardModelType const& rewardModel, std::vector const* exitRateVector); + static std::vector computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, RewardModelType const& rewardModel, std::vector const* exitRateVector); template - static std::vector computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, std::vector const& stateRewardVector, std::vector const* exitRateVector); + static std::vector computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, std::vector const& stateRewardVector, std::vector const* exitRateVector); template static std::vector computeReachabilityTimes(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& targetStates, bool qualitative); @@ -119,7 +123,14 @@ namespace storm { private: template - static std::vector computeLongRunAverages(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& probabilityMatrix, std::function const& valueGetter, std::vector const* exitRateVector); + static std::vector computeLongRunAverages(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, std::function const& valueGetter, std::vector const* exitRateVector); + template + static ValueType computeLongRunAveragesForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector); + template + static ValueType computeLongRunAveragesForBsccVi(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector); + template + static ValueType computeLongRunAveragesForBsccEqSys(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const& valueGetter, std::vector const* exitRateVector); + }; } } diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index b3e8cd8c3..857a8d478 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -1375,30 +1375,38 @@ namespace storm { } #endif + template + template + ResultValueType SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, index_type const& row) const { + typename storm::storage::SparseMatrix::const_iterator it1 = this->begin(row); + typename storm::storage::SparseMatrix::const_iterator ite1 = this->end(row); + typename storm::storage::SparseMatrix::const_iterator it2 = otherMatrix.begin(row); + typename storm::storage::SparseMatrix::const_iterator ite2 = otherMatrix.end(row); + + ResultValueType result = storm::utility::zero(); + for (;it1 != ite1 && it2 != ite2; ++it1) { + if (it1->getColumn() < it2->getColumn()) { + continue; + } else { + // If the precondition of this method (i.e. that the given matrix is a submatrix + // of the current one) was fulfilled, we know now that the two elements are in + // the same column, so we can multiply and add them to the row sum vector. + STORM_LOG_ASSERT(it1->getColumn() == it2->getColumn(), "The given matrix is not a submatrix of this one."); + result += it2->getValue() * OtherValueType(it1->getValue()); + ++it2; + } + } + return result; + } + template template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const { - std::vector result(rowCount, storm::utility::zero()); - - // Iterate over all elements of the current matrix and either continue with the next element in case the - // given matrix does not have a non-zero element at this column position, or multiply the two entries and - // add the result to the corresponding position in the vector. + std::vector result; + result.reserve(rowCount); for (index_type row = 0; row < rowCount && row < otherMatrix.getRowCount(); ++row) { - typename storm::storage::SparseMatrix::const_iterator it2 = otherMatrix.begin(row); - typename storm::storage::SparseMatrix::const_iterator ite2 = otherMatrix.end(row); - for (const_iterator it1 = this->begin(row), ite1 = this->end(row); it1 != ite1 && it2 != ite2; ++it1) { - if (it1->getColumn() < it2->getColumn()) { - continue; - } else { - // If the precondition of this method (i.e. that the given matrix is a submatrix - // of the current one) was fulfilled, we know now that the two elements are in - // the same column, so we can multiply and add them to the row sum vector. - result[row] += it2->getValue() * OtherValueType(it1->getValue()); - ++it2; - } - } + result.push_back(getPointwiseProductRowSum(otherMatrix, row)); } - return result; } diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index 125b3bd39..880c1ee36 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -819,6 +819,19 @@ namespace storm { */ std::pair, std::vector> getJacobiDecomposition() const; + /*! + * Performs a pointwise multiplication of the entries in the given row of this matrix and the entries of + * the given row of the other matrix and returns the sum. + * + * @param otherMatrix A reference to the matrix with which to perform the pointwise multiplication. This + * matrix must be a submatrix of the current matrix in the sense that it may not have entries at indices + * where there is no entry in the current matrix. + * @return the sum of the product of the entries in the given row. + */ + template + ResultValueType getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, index_type const& row) const; + + /*! * Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a vector * containing the sum of the entries in each row of the resulting matrix.