From 0ffbda5aff7eb50d7931d286b2dd7271c984fe1e Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 25 Dec 2015 17:36:37 +0100 Subject: [PATCH] initial draft of long-run rewards for parametric models Former-commit-id: 991512a57d6370ffa256a80e9d472e00c82e6eef --- examples/pdtmc/tiny/tiny.pm | 9 +- src/builder/ExplicitPrismModelBuilder.cpp | 3 - src/modelchecker/AbstractModelChecker.cpp | 2 +- .../SparseDtmcEliminationModelChecker.cpp | 111 +++++++++++++----- .../SparseDtmcEliminationModelChecker.h | 2 +- 5 files changed, 91 insertions(+), 36 deletions(-) diff --git a/examples/pdtmc/tiny/tiny.pm b/examples/pdtmc/tiny/tiny.pm index fca4f4de3..e644a10fc 100644 --- a/examples/pdtmc/tiny/tiny.pm +++ b/examples/pdtmc/tiny/tiny.pm @@ -1,11 +1,16 @@ dtmc module tiny - s : [0 .. 2] init 0; + s : [0 .. 3] init 0; - [] s = 0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [] s = 0 -> 1/3 : (s'=1) + 1/3 : (s'=2) + 1/3 : (s'=3); [] s = 1 -> 1 : (s'=2); [] s = 2 -> 1/2 : (s'=2) + 1/2 : (s'=1); + [] s = 3 -> 1 : (s'=3); endmodule +rewards + s=1 : 10; + s=3 : 5; +endrewards diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp index cb02d0d24..abf5638b3 100644 --- a/src/builder/ExplicitPrismModelBuilder.cpp +++ b/src/builder/ExplicitPrismModelBuilder.cpp @@ -1016,9 +1016,6 @@ namespace storm { for (auto const& bitVectorIndexPair : internalStateInformation.stateStorage) { stateInformation.get().valuations[bitVectorIndexPair.second] = unpackStateIntoValuation(bitVectorIndexPair.first, variableInformation); } - for (auto const& el : stateInformation.get().valuations) { - std::cout << "state: " << el << std::endl; - } } return modelComponents; diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp index 0b0d63169..2b840d83c 100644 --- a/src/modelchecker/AbstractModelChecker.cpp +++ b/src/modelchecker/AbstractModelChecker.cpp @@ -72,7 +72,7 @@ namespace storm { return this->computeInstantaneousRewards(rewardPathFormula.asInstantaneousRewardFormula(), rewardModelName, qualitative, optimalityType); } else if (rewardPathFormula.isReachabilityRewardFormula()) { return this->computeReachabilityRewards(rewardPathFormula.asReachabilityRewardFormula(), rewardModelName, qualitative, optimalityType); - } else if (rewardPathFormula.isLongRunAverageOperatorFormula()) { + } else if (rewardPathFormula.isLongRunAverageRewardFormula()) { return this->computeLongRunAverageRewards(rewardPathFormula.asLongRunAverageRewardFormula(), rewardModelName, qualitative, optimalityType); } STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << rewardPathFormula << "' is invalid."); diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index 0acec057f..6ec1e6d86 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -164,11 +164,9 @@ namespace storm { } if (furtherComputationNeeded) { - storm::storage::BitVector allStates(transitionMatrix.getRowCount(), true); - if (computeResultsForInitialStatesOnly) { // Determine the set of states that is reachable from the initial state without jumping over a target state. - storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStates, allStates, allStates, true); + storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStates, storm::storage::BitVector(numberOfStates, true), storm::storage::BitVector(numberOfStates, false)); // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state). maybeStates &= reachableStates; @@ -176,7 +174,7 @@ namespace storm { std::vector stateValues(maybeStates.size(), storm::utility::zero()); storm::utility::vector::setVectorValues(stateValues, psiStates, storm::utility::one()); - result = computeLongRunValues(transitionMatrix, this->getModel().getBackwardTransitions(), initialStates, maybeStates, computeResultsForInitialStatesOnly, stateValues); + result = computeLongRunValues(transitionMatrix, backwardTransitions, initialStates, maybeStates, computeResultsForInitialStatesOnly, stateValues); } // Construct check result based on whether we have computed values for all states or just the initial states. @@ -191,26 +189,81 @@ namespace storm { template std::unique_ptr SparseDtmcEliminationModelChecker::computeLongRunAverageRewards(storm::logic::LongRunAverageRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { - // FIXME - return nullptr; + // Do some sanity checks to establish some required properties. + RewardModelType const& rewardModel = this->getModel().getRewardModel(rewardModelName ? rewardModelName.get() : ""); + STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::IllegalArgumentException, "Input model does not have a reward model."); + + storm::storage::BitVector const& initialStates = this->getModel().getInitialStates(); + STORM_LOG_THROW(initialStates.getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); + STORM_LOG_THROW(this->computeResultsForInitialStatesOnly, storm::exceptions::IllegalArgumentException, "Cannot compute long-run probabilities for all states."); + + storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); + uint_fast64_t numberOfStates = transitionMatrix.getRowCount(); + + // Get the state-reward values from the reward model. + std::vector stateRewardValues = rewardModel.getTotalRewardVector(this->getModel().getTransitionMatrix()); + + storm::storage::BitVector maybeStates(stateRewardValues.size()); + uint_fast64_t index = 0; + for (auto const& value : stateRewardValues) { + if (value != storm::utility::zero()) { + maybeStates.set(index, true); + } + ++index; + } + + storm::storage::SparseMatrix backwardTransitions = this->getModel().getBackwardTransitions(); + + storm::storage::BitVector allStates(numberOfStates, true); + maybeStates = storm::utility::graph::performProbGreater0(backwardTransitions, allStates, maybeStates); + + std::vector result(numberOfStates, storm::utility::zero()); + + // Determine whether we need to perform some further computation. + bool furtherComputationNeeded = true; + if (computeResultsForInitialStatesOnly && initialStates.isDisjointFrom(maybeStates)) { + furtherComputationNeeded = false; + } + + if (furtherComputationNeeded) { + if (computeResultsForInitialStatesOnly) { + // Determine the set of states that is reachable from the initial state without jumping over a target state. + storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStates, storm::storage::BitVector(numberOfStates, true), storm::storage::BitVector(numberOfStates, false)); + + // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state). + maybeStates &= reachableStates; + } + + result = computeLongRunValues(transitionMatrix, backwardTransitions, initialStates, maybeStates, computeResultsForInitialStatesOnly, stateRewardValues); + } + + // Construct check result based on whether we have computed values for all states or just the initial states. + std::unique_ptr checkResult(new ExplicitQuantitativeCheckResult(result)); + if (computeResultsForInitialStatesOnly) { + // If we computed the results for the initial states only, we need to filter the result to only + // communicate these results. + checkResult->filter(ExplicitQualitativeCheckResult(initialStates)); + } + return checkResult; } template - std::vector::ValueType> SparseDtmcEliminationModelChecker::computeLongRunValues(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& statesWithProbabilityGreater0, bool computeResultsForInitialStatesOnly, std::vector& stateValues) { + std::vector::ValueType> SparseDtmcEliminationModelChecker::computeLongRunValues(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& maybeStates, bool computeResultsForInitialStatesOnly, std::vector& stateValues) { std::chrono::high_resolution_clock::time_point totalTimeStart = std::chrono::high_resolution_clock::now(); // Start by decomposing the DTMC into its BSCCs. - // FIXME: time this as well. + std::chrono::high_resolution_clock::time_point sccDecompositionStart = std::chrono::high_resolution_clock::now(); storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowCount(), true), false, true); - + auto sccDecompositionEnd = std::chrono::high_resolution_clock::now(); + std::chrono::high_resolution_clock::time_point conversionStart = std::chrono::high_resolution_clock::now(); // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(transitionMatrix); - flexibleMatrix.filter(statesWithProbabilityGreater0, statesWithProbabilityGreater0); + flexibleMatrix.filter(maybeStates, maybeStates); FlexibleSparseMatrix flexibleBackwardTransitions = getFlexibleSparseMatrix(backwardTransitions); - flexibleBackwardTransitions.filter(statesWithProbabilityGreater0, statesWithProbabilityGreater0); + flexibleBackwardTransitions.filter(maybeStates, maybeStates); auto conversionEnd = std::chrono::high_resolution_clock::now(); std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now(); @@ -223,7 +276,7 @@ namespace storm { } uint_fast64_t numberOfStates = transitionMatrix.getRowCount(); - storm::storage::BitVector statesInBsccs(numberOfStates); + storm::storage::BitVector regularStatesInBsccs(numberOfStates); storm::storage::BitVector relevantBsccs(bsccDecomposition.size()); storm::storage::BitVector bsccRepresentativesAsBitVector(numberOfStates); std::vector bsccRepresentatives; @@ -231,28 +284,24 @@ namespace storm { for (auto const& bscc : bsccDecomposition) { // Since all states in an SCC can reach all other states, we only need to check whether an arbitrary // state is a maybe state. - if (statesWithProbabilityGreater0.get(*bscc.cbegin())) { + if (maybeStates.get(*bscc.cbegin())) { relevantBsccs.set(currentIndex); bsccRepresentatives.push_back(*bscc.cbegin()); bsccRepresentativesAsBitVector.set(*bscc.cbegin(), true); for (auto const& state : bscc) { - statesInBsccs.set(state, true); + regularStatesInBsccs.set(state, true); } } ++currentIndex; } - statesInBsccs &= ~bsccRepresentativesAsBitVector; + regularStatesInBsccs &= ~bsccRepresentativesAsBitVector; // Compute the average time to stay in each state for all states in BSCCs. std::vector averageTimeInStates(stateValues.size(), storm::utility::one()); - for (auto const& el : stateValues) { - std::cout << el << std::endl; - } - // First, we eliminate all states in BSCCs (except for the representative states). { - std::unique_ptr priorityQueue = createStatePriorityQueue(distanceBasedPriorities, flexibleMatrix, flexibleBackwardTransitions, stateValues, statesInBsccs); + std::unique_ptr priorityQueue = createStatePriorityQueue(distanceBasedPriorities, flexibleMatrix, flexibleBackwardTransitions, stateValues, regularStatesInBsccs); ValueUpdateCallback valueUpdateCallback = [&stateValues,&averageTimeInStates] (storm::storage::sparse::state_type const& state, ValueType const& loopProbability) { stateValues[state] = storm::utility::simplify(loopProbability * stateValues[state]); @@ -269,17 +318,12 @@ namespace storm { }); boost::optional predecessorFilterCallback = boost::none; -// PredecessorFilterCallback([&statesInBsccs,&bsccRepresentativesAsBitVector] (storm::storage::sparse::state_type const& state) { return statesInBsccs.get(state) || bsccRepresentativesAsBitVector.get(state); } ); while (priorityQueue->hasNextState()) { storm::storage::sparse::state_type state = priorityQueue->popNextState(); eliminateState(state, flexibleMatrix, flexibleBackwardTransitions, valueUpdateCallback, predecessorCallback, priorityUpdateCallback, predecessorFilterCallback, true); STORM_LOG_ASSERT(checkConsistent(flexibleMatrix, flexibleBackwardTransitions), "The forward and backward transition matrices became inconsistent."); } - flexibleMatrix.print(); - for (auto const& el : stateValues) { - std::cout << el << std::endl; - } } // Now, we set the values of all states in BSCCs to that of the representative value (and clear the @@ -318,21 +362,29 @@ namespace storm { } // If there are states remaining that are not in BSCCs, we need to eliminate them now. - storm::storage::BitVector remainingStates = statesWithProbabilityGreater0 & ~statesInBsccs; + storm::storage::BitVector remainingStates = maybeStates & ~regularStatesInBsccs; - // Reset the values of all remaining non-BSCC states (they might have been erroneously changed by the - // previous state elimination. + // Set the value initial value of all states not in a BSCC to zero, because a) any previous value would + // incorrectly influence the result and b) the value have been erroneously changed for the predecessors of + // BSCCs by the previous state elimination. for (auto state : remainingStates) { if (!bsccRepresentativesAsBitVector.get(state)) { stateValues[state] = storm::utility::zero(); } } - performOrdinaryStateElimination(flexibleMatrix, flexibleBackwardTransitions, remainingStates, initialStates, computeResultsForInitialStatesOnly, stateValues, distanceBasedPriorities); + + // We only need to eliminate the remaining states if there was some BSCC that has a non-zero value, i.e. + // that consists of maybe states. + if (!relevantBsccs.empty()) { + performOrdinaryStateElimination(flexibleMatrix, flexibleBackwardTransitions, remainingStates, initialStates, computeResultsForInitialStatesOnly, stateValues, distanceBasedPriorities); + } std::chrono::high_resolution_clock::time_point modelCheckingEnd = std::chrono::high_resolution_clock::now(); std::chrono::high_resolution_clock::time_point totalTimeEnd = std::chrono::high_resolution_clock::now(); if (storm::settings::generalSettings().isShowStatisticsSet()) { + std::chrono::high_resolution_clock::duration sccDecompositionTime = sccDecompositionEnd - sccDecompositionStart; + std::chrono::milliseconds sccDecompositionTimeInMilliseconds = std::chrono::duration_cast(sccDecompositionTime); std::chrono::high_resolution_clock::duration conversionTime = conversionEnd - conversionStart; std::chrono::milliseconds conversionTimeInMilliseconds = std::chrono::duration_cast(conversionTime); std::chrono::high_resolution_clock::duration modelCheckingTime = modelCheckingEnd - modelCheckingStart; @@ -342,6 +394,7 @@ namespace storm { STORM_PRINT_AND_LOG(std::endl); STORM_PRINT_AND_LOG("Time breakdown:" << std::endl); + STORM_PRINT_AND_LOG(" * time for SCC decomposition: " << sccDecompositionTimeInMilliseconds.count() << "ms" << std::endl); STORM_PRINT_AND_LOG(" * time for conversion: " << conversionTimeInMilliseconds.count() << "ms" << std::endl); STORM_PRINT_AND_LOG(" * time for checking: " << modelCheckingTimeInMilliseconds.count() << "ms" << std::endl); STORM_PRINT_AND_LOG("------------------------------------------" << std::endl); diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h index d8eda980c..692a6be57 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h @@ -113,7 +113,7 @@ namespace storm { PenaltyFunctionType penaltyFunction; }; - static std::vector computeLongRunValues(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& statesWithProbabilityGreater0, bool computeResultsForInitialStatesOnly, std::vector& stateValues); + static std::vector computeLongRunValues(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& maybeStates, bool computeResultsForInitialStatesOnly, std::vector& stateValues); static std::vector computeReachabilityValues(storm::storage::SparseMatrix const& transitionMatrix, std::vector& values, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, bool computeResultsForInitialStatesOnly, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& oneStepProbabilitiesToTarget);