From 1860502a3afb4be68256155fff88219051f5cb90 Mon Sep 17 00:00:00 2001 From: TimQu Date: Mon, 19 Oct 2015 19:56:28 +0200 Subject: [PATCH] Deterministic states with only constant outgoing transitions are now eliminated Former-commit-id: be5bf4f7ccffe941131e5d7cf44d89d564d39423 --- .../SparseDtmcEliminationModelChecker.cpp | 335 +----------------- .../SparseDtmcEliminationModelChecker.h | 39 +- .../region/ApproximationModel.cpp | 2 +- src/modelchecker/region/SamplingModel.cpp | 2 +- .../region/SparseDtmcRegionModelChecker.cpp | 20 +- .../region/SparseDtmcRegionModelChecker.h | 4 - .../region/SparseMdpRegionModelChecker.cpp | 116 ++++-- src/storage/FlexibleSparseMatrix.cpp | 326 +++++++++++++++++ src/storage/FlexibleSparseMatrix.h | 51 +++ .../SparseMdpRegionModelCheckerTest.cpp | 2 +- 10 files changed, 505 insertions(+), 392 deletions(-) create mode 100644 src/storage/FlexibleSparseMatrix.cpp create mode 100644 src/storage/FlexibleSparseMatrix.h diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index aefc6f023..d7c1d06ed 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -268,18 +268,18 @@ namespace storm { STORM_LOG_INFO("Eliminating " << states.size() << " states using the state elimination technique." << std::endl); boost::optional> missingStateRewards; std::chrono::high_resolution_clock::time_point conversionStart = std::chrono::high_resolution_clock::now(); - FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(submatrix); - FlexibleSparseMatrix flexibleBackwardTransitions = getFlexibleSparseMatrix(submatrixTransposed, true); + storm::storage::FlexibleSparseMatrix flexibleMatrix(submatrix); + storm::storage::FlexibleSparseMatrix flexibleBackwardTransitions(submatrixTransposed, true); std::chrono::high_resolution_clock::time_point conversionEnd = std::chrono::high_resolution_clock::now(); std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now(); for (auto const& state : states) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); } STORM_LOG_INFO("Eliminated " << states.size() << " states." << std::endl); // Eliminate the transitions going into the initial state (if there are any). if (!flexibleBackwardTransitions.getRow(*newInitialStates.begin()).empty()) { - eliminateState(flexibleMatrix, oneStepProbabilities, *newInitialStates.begin(), flexibleBackwardTransitions, missingStateRewards, false); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, *newInitialStates.begin(), flexibleBackwardTransitions, missingStateRewards, false); } // Now we need to basically eliminate all chains of not-psi states after phi states and chains of not-phi @@ -310,11 +310,11 @@ namespace storm { for (auto const& element : currentRow) { // If any of the successors is a phi state, we eliminate it (wrt. all its phi predecessors). if (!psiStates.get(element.getColumn())) { - typename FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); + typename storm::storage::FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); // Eliminate the successor only if there possibly is a psi state reachable through it. if (successorRow.size() > 1 || (!successorRow.empty() && successorRow.front().getColumn() != element.getColumn())) { STORM_LOG_TRACE("Found non-psi successor " << element.getColumn() << " that needs to be eliminated."); - eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, phiStates); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, phiStates); hasNonPsiSuccessor = true; } } @@ -340,10 +340,10 @@ namespace storm { for (auto const& element : currentRow) { // If any of the successors is a psi state, we eliminate it (wrt. all its psi predecessors). if (!phiStates.get(element.getColumn())) { - typename FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); + typename storm::storage::FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); if (successorRow.size() > 1 || (!successorRow.empty() && successorRow.front().getColumn() != element.getColumn())) { STORM_LOG_TRACE("Found non-phi successor " << element.getColumn() << " that needs to be eliminated."); - eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, psiStates); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, psiStates); hasNonPhiSuccessor = true; } } @@ -417,8 +417,8 @@ namespace storm { 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); - FlexibleSparseMatrix flexibleBackwardTransitions = getFlexibleSparseMatrix(backwardTransitions, true); + storm::storage::FlexibleSparseMatrix flexibleMatrix(transitionMatrix); + storm::storage::FlexibleSparseMatrix flexibleBackwardTransitions(backwardTransitions, true); auto conversionEnd = std::chrono::high_resolution_clock::now(); std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now(); @@ -437,7 +437,7 @@ namespace storm { STORM_LOG_DEBUG("Eliminating " << states.size() << " states using the state elimination technique." << std::endl); for (auto const& state : states) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); } STORM_LOG_DEBUG("Eliminated " << states.size() << " states." << std::endl); } else if (storm::settings::sparseDtmcEliminationModelCheckerSettings().getEliminationMethod() == storm::settings::modules::SparseDtmcEliminationModelCheckerSettings::EliminationMethod::Hybrid) { @@ -450,7 +450,7 @@ namespace storm { STORM_LOG_DEBUG("Eliminating " << entryStateQueue.size() << " entry states as a last step."); if (storm::settings::sparseDtmcEliminationModelCheckerSettings().isEliminateEntryStatesLastSet()) { for (auto const& state : entryStateQueue) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); } } STORM_LOG_DEBUG("Eliminated " << subsystem.size() << " states." << std::endl); @@ -461,7 +461,7 @@ namespace storm { // If we are computing probabilities, then we can simply call the state elimination procedure. It // will scale the transition row of the initial state with 1/(1-loopProbability). STORM_LOG_INFO("Eliminating initial state " << *initialStates.begin() << "." << std::endl); - eliminateState(flexibleMatrix, oneStepProbabilities, *initialStates.begin(), flexibleBackwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleMatrix, oneStepProbabilities, *initialStates.begin(), flexibleBackwardTransitions, stateRewards); } else { // If we are computing rewards, we cannot call the state elimination procedure for technical reasons. // Instead, we need to get rid of a potential loop in this state explicitly. @@ -563,7 +563,7 @@ namespace storm { } template - uint_fast64_t SparseDtmcEliminationModelChecker::treatScc(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, boost::optional>& stateRewards, boost::optional> const& statePriorities) { + uint_fast64_t SparseDtmcEliminationModelChecker::treatScc(storm::storage::FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, storm::storage::FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, boost::optional>& stateRewards, boost::optional> const& statePriorities) { uint_fast64_t maximalDepth = level; // If the SCCs are large enough, we try to split them further. @@ -595,7 +595,7 @@ namespace storm { STORM_LOG_TRACE("Eliminating " << trivialSccs.size() << " trivial SCCs."); for (auto const& stateIndexPair : trivialSccs) { - eliminateState(matrix, oneStepProbabilities, stateIndexPair.first, backwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(matrix, oneStepProbabilities, stateIndexPair.first, backwardTransitions, stateRewards); remainingSccs.set(stateIndexPair.second, false); } STORM_LOG_TRACE("Eliminated all trivial SCCs."); @@ -638,7 +638,7 @@ namespace storm { // Eliminate the remaining states that do not have a self-loop (in the current, i.e. modified) // transition probability matrix. for (auto const& state : states) { - eliminateState(matrix, oneStepProbabilities, state, backwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(matrix, oneStepProbabilities, state, backwardTransitions, stateRewards); } STORM_LOG_TRACE("Eliminated all states of SCC."); @@ -648,7 +648,7 @@ namespace storm { if (eliminateEntryStates) { STORM_LOG_TRACE("Finally, eliminating/adding entry states."); for (auto state : entryStates) { - eliminateState(matrix, oneStepProbabilities, state, backwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(matrix, oneStepProbabilities, state, backwardTransitions, stateRewards); } STORM_LOG_TRACE("Eliminated/added entry states."); } else { @@ -660,307 +660,6 @@ namespace storm { return maximalDepth; } - template - void SparseDtmcEliminationModelChecker::eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix& backwardTransitions, boost::optional>& stateRewards, bool removeForwardTransitions, bool constrained, storm::storage::BitVector const& predecessorConstraint) { - - bool hasSelfLoop = false; - ValueType loopProbability = storm::utility::zero(); - - // Start by finding loop probability. - typename FlexibleSparseMatrix::row_type& currentStateSuccessors = matrix.getRow(state); - for (auto entryIt = currentStateSuccessors.begin(), entryIte = currentStateSuccessors.end(); entryIt != entryIte; ++entryIt) { - if (entryIt->getColumn() >= state) { - if (entryIt->getColumn() == state) { - loopProbability = entryIt->getValue(); - hasSelfLoop = true; - - // If we do not clear the forward transitions completely, we need to remove the self-loop, - // because we scale all the other outgoing transitions with it anyway.. - if (!removeForwardTransitions) { - currentStateSuccessors.erase(entryIt); - } - } - break; - } - } - - // Scale all entries in this row with (1 / (1 - loopProbability)) only in case there was a self-loop. - std::size_t scaledSuccessors = 0; - if (hasSelfLoop) { - STORM_LOG_ASSERT(loopProbability != storm::utility::one(), "Must not eliminate state with probability 1 self-loop."); - loopProbability = storm::utility::one() / (storm::utility::one() - loopProbability); - storm::utility::simplify(loopProbability); - for (auto& entry : matrix.getRow(state)) { - // Only scale the non-diagonal entries. - if (entry.getColumn() != state) { - ++scaledSuccessors; - entry.setValue(storm::utility::simplify(entry.getValue() * loopProbability)); - } - } - if (!stateRewards) { - oneStepProbabilities[state] = oneStepProbabilities[state] * loopProbability; - } - } - - STORM_LOG_TRACE((hasSelfLoop ? "State has self-loop." : "State does not have a self-loop.")); - - // Now connect the predecessors of the state being eliminated with its successors. - typename FlexibleSparseMatrix::row_type& currentStatePredecessors = backwardTransitions.getRow(state); - std::size_t predecessorForwardTransitionCount = 0; - - // In case we have a constrained elimination, we need to keep track of the new predecessors. - typename FlexibleSparseMatrix::row_type newCurrentStatePredecessors; - - // Now go through the predecessors and eliminate the ones (satisfying the constraint if given). - for (auto const& predecessorEntry : currentStatePredecessors) { - uint_fast64_t predecessor = predecessorEntry.getColumn(); - - // Skip the state itself as one of its predecessors. - if (predecessor == state) { - assert(hasSelfLoop); - continue; - } - - // Skip the state if the elimination is constrained, but the predecessor is not in the constraint. - if (constrained && !predecessorConstraint.get(predecessor)) { - newCurrentStatePredecessors.emplace_back(predecessor, storm::utility::one()); - STORM_LOG_TRACE("Not eliminating predecessor " << predecessor << ", because it does not fit the filter."); - continue; - } - STORM_LOG_TRACE("Eliminating predecessor " << predecessor << "."); - - // First, find the probability with which the predecessor can move to the current state, because - // the other probabilities need to be scaled with this factor. - typename FlexibleSparseMatrix::row_type& predecessorForwardTransitions = matrix.getRow(predecessor); - predecessorForwardTransitionCount += predecessorForwardTransitions.size(); - typename FlexibleSparseMatrix::row_type::iterator multiplyElement = std::find_if(predecessorForwardTransitions.begin(), predecessorForwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == state; }); - - // Make sure we have found the probability and set it to zero. - STORM_LOG_THROW(multiplyElement != predecessorForwardTransitions.end(), storm::exceptions::InvalidStateException, "No probability for successor found."); - ValueType multiplyFactor = multiplyElement->getValue(); - multiplyElement->setValue(storm::utility::zero()); - - // At this point, we need to update the (forward) transitions of the predecessor. - typename FlexibleSparseMatrix::row_type::iterator first1 = predecessorForwardTransitions.begin(); - typename FlexibleSparseMatrix::row_type::iterator last1 = predecessorForwardTransitions.end(); - typename FlexibleSparseMatrix::row_type::iterator first2 = currentStateSuccessors.begin(); - typename FlexibleSparseMatrix::row_type::iterator last2 = currentStateSuccessors.end(); - - typename FlexibleSparseMatrix::row_type newSuccessors; - newSuccessors.reserve((last1 - first1) + (last2 - first2)); - std::insert_iterator result(newSuccessors, newSuccessors.end()); - - // Now we merge the two successor lists. (Code taken from std::set_union and modified to suit our needs). - for (; first1 != last1; ++result) { - // Skip the transitions to the state that is currently being eliminated. - if (first1->getColumn() == state || (first2 != last2 && first2->getColumn() == state)) { - if (first1->getColumn() == state) { - ++first1; - } - if (first2 != last2 && first2->getColumn() == state) { - ++first2; - } - continue; - } - - if (first2 == last2) { - std::copy_if(first1, last1, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; } ); - break; - } - if (first2->getColumn() < first1->getColumn()) { - *result = storm::utility::simplify(std::move(*first2 * multiplyFactor)); - ++first2; - } else if (first1->getColumn() < first2->getColumn()) { - *result = *first1; - ++first1; - } else { - *result = storm::storage::MatrixEntry(first1->getColumn(), storm::utility::simplify(first1->getValue() + storm::utility::simplify(multiplyFactor * first2->getValue()))); - ++first1; - ++first2; - } - } - for (; first2 != last2; ++first2) { - if (first2->getColumn() != state) { - *result = storm::utility::simplify(std::move(*first2 * multiplyFactor)); - } - } - - // Now move the new transitions in place. - predecessorForwardTransitions = std::move(newSuccessors); - - if (!stateRewards) { - // Add the probabilities to go to a target state in just one step if we have to compute probabilities. - oneStepProbabilities[predecessor] += storm::utility::simplify(multiplyFactor * oneStepProbabilities[state]); - STORM_LOG_TRACE("Fixed new next-state probabilities of predecessor states."); - } else { - // If we are computing rewards, we basically scale the state reward of the state to eliminate and - // add the result to the state reward of the predecessor. - if (hasSelfLoop) { - stateRewards.get()[predecessor] += storm::utility::simplify(multiplyFactor * loopProbability * stateRewards.get()[state]); - } else { - stateRewards.get()[predecessor] += storm::utility::simplify(multiplyFactor * stateRewards.get()[state]); - } - } - } - - // Finally, we need to add the predecessor to the set of predecessors of every successor. - for (auto const& successorEntry : currentStateSuccessors) { - typename FlexibleSparseMatrix::row_type& successorBackwardTransitions = backwardTransitions.getRow(successorEntry.getColumn()); - - // Delete the current state as a predecessor of the successor state only if we are going to remove the - // current state's forward transitions. - if (removeForwardTransitions) { - typename FlexibleSparseMatrix::row_type::iterator elimIt = std::find_if(successorBackwardTransitions.begin(), successorBackwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == state; }); - STORM_LOG_ASSERT(elimIt != successorBackwardTransitions.end(), "Expected a proper backward transition, but found none."); - successorBackwardTransitions.erase(elimIt); - } - - typename FlexibleSparseMatrix::row_type::iterator first1 = successorBackwardTransitions.begin(); - typename FlexibleSparseMatrix::row_type::iterator last1 = successorBackwardTransitions.end(); - typename FlexibleSparseMatrix::row_type::iterator first2 = currentStatePredecessors.begin(); - typename FlexibleSparseMatrix::row_type::iterator last2 = currentStatePredecessors.end(); - - typename FlexibleSparseMatrix::row_type newPredecessors; - newPredecessors.reserve((last1 - first1) + (last2 - first2)); - std::insert_iterator result(newPredecessors, newPredecessors.end()); - - if (!constrained) { - for (; first1 != last1; ++result) { - if (first2 == last2) { - std::copy(first1, last1, result); - break; - } - if (first2->getColumn() < first1->getColumn()) { - if (first2->getColumn() != state) { - *result = *first2; - } - ++first2; - } else { - *result = *first1; - if (first1->getColumn() == first2->getColumn()) { - ++first2; - } - ++first1; - } - } - std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; }); - } else { - // If the elimination is constrained, we need to be more selective when we set the new predecessors - // of the successor state. - for (; first1 != last1; ++result) { - if (first2 == last2) { - std::copy(first1, last1, result); - break; - } - if (first2->getColumn() < first1->getColumn()) { - if (first2->getColumn() != state) { - *result = *first2; - } - ++first2; - } else { - *result = *first1; - if (first1->getColumn() == first2->getColumn()) { - ++first2; - } - ++first1; - } - } - std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state && (!constrained || predecessorConstraint.get(a.getColumn())); }); - } - - // Now move the new predecessors in place. - successorBackwardTransitions = std::move(newPredecessors); - } - STORM_LOG_TRACE("Fixed predecessor lists of successor states."); - - if (removeForwardTransitions) { - // Clear the eliminated row to reduce memory consumption. - currentStateSuccessors.clear(); - currentStateSuccessors.shrink_to_fit(); - } - if (!constrained) { - currentStatePredecessors.clear(); - currentStatePredecessors.shrink_to_fit(); - } else { - currentStatePredecessors = std::move(newCurrentStatePredecessors); - } - } - - template - SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::FlexibleSparseMatrix(index_type rows) : data(rows) { - // Intentionally left empty. - } - - template - void SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::reserveInRow(index_type row, index_type numberOfElements) { - this->data[row].reserve(numberOfElements); - } - - template - typename SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::row_type& SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::getRow(index_type index) { - return this->data[index]; - } - - template - typename SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::row_type const& SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::getRow(index_type index) const { - return this->data[index]; - } - - template - typename SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::index_type SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::getNumberOfRows() const { - return this->data.size(); - } - - template - bool SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::hasSelfLoop(storm::storage::sparse::state_type state) { - for (auto const& entry : this->getRow(state)) { - if (entry.getColumn() < state) { - continue; - } else if (entry.getColumn() > state) { - return false; - } else if (entry.getColumn() == state) { - return true; - } - } - return false; - } - - template - void SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::print() const { - for (uint_fast64_t index = 0; index < this->data.size(); ++index) { - std::cout << index << " - "; - for (auto const& element : this->getRow(index)) { - std::cout << "(" << element.getColumn() << ", " << element.getValue() << ") "; - } - std::cout << std::endl; - } - } - - template - typename SparseDtmcEliminationModelChecker::FlexibleSparseMatrix SparseDtmcEliminationModelChecker::getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne) { - FlexibleSparseMatrix flexibleMatrix(matrix.getRowCount()); - - for (typename FlexibleSparseMatrix::index_type rowIndex = 0; rowIndex < matrix.getRowCount(); ++rowIndex) { - typename storm::storage::SparseMatrix::const_rows row = matrix.getRow(rowIndex); - flexibleMatrix.reserveInRow(rowIndex, row.getNumberOfEntries()); - - for (auto const& element : row) { - // If the probability is zero, we skip this entry. - if (storm::utility::isZero(element.getValue())) { - continue; - } - - if (setAllValuesToOne) { - flexibleMatrix.getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one()); - } else { - flexibleMatrix.getRow(rowIndex).emplace_back(element); - } - } - } - - return flexibleMatrix; - } - template class SparseDtmcEliminationModelChecker>; #ifdef STORM_HAVE_CARL diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h index d2ad1570f..1d712b4db 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h @@ -2,6 +2,7 @@ #define STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCELIMINATIONMODELCHECKER_H_ #include "src/storage/sparse/StateType.h" +#include "src/storage/FlexibleSparseMatrix.h" #include "src/models/sparse/Dtmc.h" #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "src/utility/constants.h" @@ -35,46 +36,10 @@ namespace storm { virtual std::unique_ptr computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; private: - class FlexibleSparseMatrix { - public: - typedef uint_fast64_t index_type; - typedef ValueType value_type; - typedef std::vector> row_type; - typedef typename row_type::iterator iterator; - typedef typename row_type::const_iterator const_iterator; - - FlexibleSparseMatrix() = default; - FlexibleSparseMatrix(index_type rows); - - void reserveInRow(index_type row, index_type numberOfElements); - - row_type& getRow(index_type); - row_type const& getRow(index_type) const; - - index_type getNumberOfRows() const; - - void print() const; - - /*! - * Checks whether the given state has a self-loop with an arbitrary probability in the given probability matrix. - * - * @param state The state for which to check whether it possesses a self-loop. - * @param matrix The matrix in which to look for the loop. - * @return True iff the given state has a self-loop with an arbitrary probability in the given probability matrix. - */ - bool hasSelfLoop(storm::storage::sparse::state_type state); - - private: - std::vector data; - }; ValueType computeReachabilityValue(storm::storage::SparseMatrix const& transitionMatrix, std::vector& oneStepProbabilities, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional>& stateRewards, boost::optional> const& statePriorities = {}); - uint_fast64_t treatScc(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, boost::optional>& stateRewards, boost::optional> const& statePriorities = {}); - - static FlexibleSparseMatrix getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false); - - void eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix& backwardTransitions, boost::optional>& stateRewards, bool removeForwardTransitions = true, bool constrained = false, storm::storage::BitVector const& predecessorConstraint = storm::storage::BitVector()); + uint_fast64_t treatScc(storm::storage::FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, storm::storage::FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, boost::optional>& stateRewards, boost::optional> const& statePriorities = {}); std::vector getStatePriorities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& transitionMatrixTransposed, storm::storage::BitVector const& initialStates, std::vector const& oneStepProbabilities); }; diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index f6265a83e..c92b4d0c1 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/src/modelchecker/region/ApproximationModel.cpp @@ -166,7 +166,7 @@ namespace storm { } } if(!this->computeRewards){ - if(storm::utility::isConstant(targetProbability)){ + if(storm::utility::isConstant(storm::utility::simplify(targetProbability))){ *vectorIt = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability)); } else { auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(targetProbability, rowSubstitutions[curRow]), dummyValue)).first; diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index 79419404c..9539c2bde 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -162,7 +162,7 @@ namespace storm { } } if(!this->computeRewards){ - if(storm::utility::isConstant(targetProbability)){ + if(storm::utility::isConstant(storm::utility::simplify(targetProbability))){ *vectorIt = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability)); } else { typename std::unordered_map::iterator functionsIt = this->functions.insert(FunctionEntry(targetProbability, dummyValue)).first; diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp index b349a475e..a1abe1e15 100644 --- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp +++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp @@ -9,11 +9,13 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/region/RegionCheckResult.h" #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" +#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h" #include "src/models/sparse/StandardRewardModel.h" #include "src/settings/SettingsManager.h" #include "src/settings/modules/RegionSettings.h" #include "src/solver/OptimizationDirection.h" #include "src/storage/sparse/StateType.h" +#include "src/storage/FlexibleSparseMatrix.h" #include "src/utility/constants.h" #include "src/utility/graph.h" #include "src/utility/macros.h" @@ -34,8 +36,7 @@ namespace storm { template SparseDtmcRegionModelChecker::SparseDtmcRegionModelChecker(ParametricSparseModelType const& model) : - AbstractSparseRegionModelChecker(model), - eliminationModelChecker(model){ + AbstractSparseRegionModelChecker(model){ //intentionally left empty } @@ -102,8 +103,8 @@ namespace storm { storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates); // Eliminate all states with only constant outgoing transitions // Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. - auto flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix); - auto flexibleBackwardTransitions= this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix.transpose(), true); + storm::storage::FlexibleSparseMatrix flexibleTransitions(submatrix); + storm::storage::FlexibleSparseMatrix flexibleBackwardTransitions(submatrix.transpose(), true); // Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated. storm::storage::BitVector subsystem(submatrix.getRowCount(), true); //The states that we consider to eliminate @@ -140,7 +141,7 @@ namespace storm { } } if(eliminateThisState){ - this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleTransitions, oneStepProbabilities, state, state, flexibleBackwardTransitions, stateRewards); subsystem.set(state,false); } } @@ -182,8 +183,8 @@ namespace storm { matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, storm::utility::simplify(oneStepProbabilities[oldStateIndex])); } //transition to sink state - if(!storm::utility::isZero(missingProbability)){ - matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, storm::utility::simplify(missingProbability)); + if(!storm::utility::isZero(storm::utility::simplify(missingProbability))){ + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, missingProbability); } } } @@ -397,8 +398,9 @@ namespace storm { if(this->isComputeRewards()){ stateRewards = simpleModel.getUniqueRewardModel()->second.getTotalRewardVector(maybeStates.getNumberOfSetBits(), simpleModel.getTransitionMatrix(), maybeStates); } - std::vector statePriorities = this->eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities); - this->reachabilityFunction=std::make_shared(this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, simpleModel.getStates("target"), stateRewards, statePriorities)); + storm::modelchecker::SparseDtmcEliminationModelChecker eliminationModelChecker(simpleModel); + std::vector statePriorities = eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities); + this->reachabilityFunction=std::make_shared(eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, simpleModel.getStates("target"), stateRewards, statePriorities)); /* std::string funcStr = " (/ " + this->reachabilityFunction->nominator().toString(false, true) + " " + this->reachabilityFunction->denominator().toString(false, true) + diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.h b/src/modelchecker/region/SparseDtmcRegionModelChecker.h index 53bbf6d09..d82df9897 100644 --- a/src/modelchecker/region/SparseDtmcRegionModelChecker.h +++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.h @@ -6,7 +6,6 @@ #include "src/models/sparse/StandardRewardModel.h" #include "src/models/sparse/Dtmc.h" #include "src/utility/region.h" -#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h" #include "src/solver/Smt2SmtSolver.h" namespace storm { @@ -119,9 +118,6 @@ namespace storm { // The function for the reachability probability (or: reachability reward) in the initial state std::shared_ptr reachabilityFunction; - // Instance of an elimination model checker to access its functions - storm::modelchecker::SparseDtmcEliminationModelChecker> eliminationModelChecker; - // the smt solver that is used to prove properties with the help of the reachabilityFunction std::shared_ptr smtSolver; diff --git a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp index 6018d22d3..f7e284773 100644 --- a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp +++ b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp @@ -14,6 +14,7 @@ #include "src/settings/modules/RegionSettings.h" #include "src/solver/OptimizationDirection.h" #include "src/storage/sparse/StateType.h" +#include "src/storage/FlexibleSparseMatrix.h" #include "src/utility/constants.h" #include "src/utility/graph.h" #include "src/utility/macros.h" @@ -77,40 +78,113 @@ namespace storm { storm::storage::BitVector maybeStates, targetStates; preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, constantResult); - //lets count the number of states without nondet choices and with constant outgoing transitions - //TODO: Remove this - storm::storage::SparseMatrixconst& matrix=this->getModel().getTransitionMatrix(); - uint_fast64_t stateCounter=0; - for (uint_fast64_t state=0; stategetModel().getNumberOfStates();++state){ - if(matrix.getRowGroupSize(state)==1){ - bool hasConstTransitions=true; - for(auto const& entry : matrix.getRowGroup(state)){ + STORM_LOG_DEBUG("Elimination of deterministic states with constant outgoing transitions is happening now."); + // 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(this->getModel().getTransitionMatrix(), this->getModel().getInitialStates(), maybeStates, targetStates); + // 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; + // Create a vector for the probabilities to go to a target state in one step. + std::vector oneStepProbabilities = this->getModel().getTransitionMatrix().getConstrainedRowGroupSumVector(maybeStates, targetStates); + // Determine the initial state of the sub-model. + storm::storage::sparse::state_type initialState = *(this->getModel().getInitialStates() % maybeStates).begin(); + // We then build the submatrix that only has the transitions of the maybe states. + storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates); + boost::optional> noStateRewards; + // Eliminate all deterministic states with only constant outgoing transitions + // Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. + storm::storage::FlexibleSparseMatrix flexibleTransitions(submatrix); + storm::storage::FlexibleSparseMatrix flexibleBackwardTransitions(submatrix.transpose(), true); + // Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated. + storm::storage::BitVector subsystem(submatrix.getRowGroupCount(), true); + //The states that we consider to eliminate + storm::storage::BitVector considerToEliminate(submatrix.getRowGroupCount(), true); + considerToEliminate.set(initialState, false); + for (auto const& state : considerToEliminate) { + bool eliminateThisState=true; + if(submatrix.getRowGroupSize(state) == 1){ + //state is deterministic. Check if outgoing transitions are constant + for(auto const& entry : submatrix.getRowGroup(state)){ if(!storm::utility::isConstant(entry.getValue())){ - hasConstTransitions = false; + eliminateThisState=false; + break; } } - if(hasConstTransitions){ - ++stateCounter; + if(!storm::utility::isConstant(oneStepProbabilities[submatrix.getRowGroupIndices()[state]])){ + eliminateThisState=false; } + } else { + eliminateThisState = false; + } + if(eliminateThisState){ + storm::storage::FlexibleSparseMatrix::eliminateState(flexibleTransitions, oneStepProbabilities, state, submatrix.getRowGroupIndices()[state], flexibleBackwardTransitions, noStateRewards); + subsystem.set(state,false); } } - std::cout << "Found that " << stateCounter << " of " << this->getModel().getNumberOfStates() << " states could be eliminated" << std::endl; - - //TODO: Actually eliminate the states... - STORM_LOG_WARN("No simplification of the original model (like elimination of constant transitions) is happening. Will just use a copy of the original model"); + STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions."); + std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl; + + //Build the simple model + STORM_LOG_DEBUG("Building the resulting simplified model."); + //The matrix. The flexibleTransitions matrix might have empty rows where states have been eliminated. + //The new matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly. + std::vector newStateIndexMap(flexibleTransitions.getNumberOfRows(), flexibleTransitions.getNumberOfRows()); //initialize with some illegal index + storm::storage::sparse::state_type newStateIndex=0; + for(auto const& state : subsystem){ + newStateIndexMap[state]=newStateIndex; + ++newStateIndex; + } + //We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead + storm::storage::sparse::state_type numStates=newStateIndex+2; + storm::storage::sparse::state_type targetState=numStates-2; + storm::storage::sparse::state_type sinkState= numStates-1; + //We can now fill in the data. + storm::storage::SparseMatrixBuilder matrixBuilder(0, numStates, 0, true, true, numStates); + std::size_t curRow = 0; + for(auto oldRowGroup : subsystem){ + matrixBuilder.newRowGroup(curRow); + for (auto oldRow = submatrix.getRowGroupIndices()[oldRowGroup]; oldRow < submatrix.getRowGroupIndices()[oldRowGroup+1]; ++oldRow){ + ParametricType missingProbability=storm::utility::region::getNewFunction(storm::utility::one()); + //go through columns: + for(auto& entry: flexibleTransitions.getRow(oldRow)){ + STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=flexibleTransitions.getNumberOfRows(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated."); + missingProbability-=entry.getValue(); + matrixBuilder.addNextValue(curRow,newStateIndexMap[entry.getColumn()], storm::utility::simplify(entry.getValue())); + } + //transition to target state + if(!storm::utility::isZero(oneStepProbabilities[oldRow])){ + missingProbability-=oneStepProbabilities[oldRow]; + matrixBuilder.addNextValue(curRow, targetState, storm::utility::simplify(oneStepProbabilities[oldRow])); + } + //transition to sink state + if(!storm::utility::isZero(storm::utility::simplify(missingProbability))){ + matrixBuilder.addNextValue(curRow, sinkState, missingProbability); + } + ++curRow; + } + } + //add self loops on the additional states (i.e., target and sink) + matrixBuilder.newRowGroup(curRow); + matrixBuilder.addNextValue(curRow, targetState, storm::utility::one()); + ++curRow; + matrixBuilder.newRowGroup(curRow); + matrixBuilder.addNextValue(curRow, sinkState, storm::utility::one()); //Get a new labeling - storm::storage::sparse::state_type numStates = this->getModel().getNumberOfStates(); - storm::storage::BitVector sinkStates = ~(maybeStates | targetStates); storm::models::sparse::StateLabeling labeling(numStates); - storm::storage::BitVector initLabel(this->getModel().getInitialStates()); + storm::storage::BitVector initLabel(numStates, false); + initLabel.set(newStateIndexMap[initialState], true); labeling.addLabel("init", std::move(initLabel)); - labeling.addLabel("target", std::move(targetStates)); - labeling.addLabel("sink", std::move(sinkStates)); + storm::storage::BitVector targetLabel(numStates, false); + targetLabel.set(targetState, true); + labeling.addLabel("target", std::move(targetLabel)); + storm::storage::BitVector sinkLabel(numStates, false); + sinkLabel.set(sinkState, true); + labeling.addLabel("sink", std::move(sinkLabel)); + //Other ingredients std::unordered_map noRewardModels; boost::optional>> noChoiceLabeling; - simpleModel = std::make_shared(this->getModel().getTransitionMatrix(), std::move(labeling), std::move(noRewardModels), std::move(noChoiceLabeling)); + simpleModel = std::make_shared(matrixBuilder.build(), std::move(labeling), std::move(noRewardModels), std::move(noChoiceLabeling)); //Get the simplified formula std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); diff --git a/src/storage/FlexibleSparseMatrix.cpp b/src/storage/FlexibleSparseMatrix.cpp new file mode 100644 index 000000000..1b1f6eb1b --- /dev/null +++ b/src/storage/FlexibleSparseMatrix.cpp @@ -0,0 +1,326 @@ +/* + * File: FlexibleSparseMatrix.cpp + * + * Created on October 19, 2015, 5:19 PM + */ + +#include + +#include "FlexibleSparseMatrix.h" +#include "src/adapters/CarlAdapter.h" +#include "src/utility/graph.h" +#include "src/utility/vector.h" +#include "src/utility/macros.h" +#include "src/exceptions/InvalidStateException.h" + +namespace storm { + namespace storage { + + template + FlexibleSparseMatrix::FlexibleSparseMatrix(index_type rows) : data(rows) { + // Intentionally left empty. + } + + template + FlexibleSparseMatrix::FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne) : data(matrix.getRowCount()) { + + for (index_type rowIndex = 0; rowIndex < matrix.getRowCount(); ++rowIndex) { + typename storm::storage::SparseMatrix::const_rows row = matrix.getRow(rowIndex); + reserveInRow(rowIndex, row.getNumberOfEntries()); + + for (auto const& element : row) { + // If the probability is zero, we skip this entry. + if (storm::utility::isZero(element.getValue())) { + continue; + } + + if (setAllValuesToOne) { + getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one()); + } else { + getRow(rowIndex).emplace_back(element); + } + } + } + } + + template + void FlexibleSparseMatrix::reserveInRow(index_type row, index_type numberOfElements) { + this->data[row].reserve(numberOfElements); + } + + template + typename FlexibleSparseMatrix::row_type& FlexibleSparseMatrix::getRow(index_type index) { + return this->data[index]; + } + + template + typename FlexibleSparseMatrix::row_type const& FlexibleSparseMatrix::getRow(index_type index) const { + return this->data[index]; + } + + template + typename FlexibleSparseMatrix::index_type FlexibleSparseMatrix::getNumberOfRows() const { + return this->data.size(); + } + + template + bool FlexibleSparseMatrix::hasSelfLoop(storm::storage::sparse::state_type state) { + for (auto const& entry : this->getRow(state)) { + if (entry.getColumn() < state) { + continue; + } else if (entry.getColumn() > state) { + return false; + } else if (entry.getColumn() == state) { + return true; + } + } + return false; + } + + template + void FlexibleSparseMatrix::print() const { + for (uint_fast64_t index = 0; index < this->data.size(); ++index) { + std::cout << index << " - "; + for (auto const& element : this->getRow(index)) { + std::cout << "(" << element.getColumn() << ", " << element.getValue() << ") "; + } + std::cout << std::endl; + } + } + + + template + void FlexibleSparseMatrix::eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, uint_fast64_t row, FlexibleSparseMatrix& backwardTransitions, boost::optional>& stateRewards, bool removeForwardTransitions, bool constrained, storm::storage::BitVector const& predecessorConstraint) { + + bool hasSelfLoop = false; + ValueType loopProbability = storm::utility::zero(); + + // Start by finding loop probability. + row_type& currentStateSuccessors = matrix.getRow(row); + for (auto entryIt = currentStateSuccessors.begin(), entryIte = currentStateSuccessors.end(); entryIt != entryIte; ++entryIt) { + if (entryIt->getColumn() >= state) { + if (entryIt->getColumn() == state) { + loopProbability = entryIt->getValue(); + hasSelfLoop = true; + + // If we do not clear the forward transitions completely, we need to remove the self-loop, + // because we scale all the other outgoing transitions with it anyway.. + if (!removeForwardTransitions) { + currentStateSuccessors.erase(entryIt); + } + } + break; + } + } + + // Scale all entries in this row with (1 / (1 - loopProbability)) only in case there was a self-loop. + std::size_t scaledSuccessors = 0; + if (hasSelfLoop) { + STORM_LOG_ASSERT(loopProbability != storm::utility::one(), "Must not eliminate state with probability 1 self-loop."); + loopProbability = storm::utility::one() / (storm::utility::one() - loopProbability); + storm::utility::simplify(loopProbability); + for (auto& entry : matrix.getRow(row)) { + // Only scale the non-diagonal entries. + if (entry.getColumn() != state) { + ++scaledSuccessors; + entry.setValue(storm::utility::simplify(entry.getValue() * loopProbability)); + } + } + if (!stateRewards) { + oneStepProbabilities[row] = oneStepProbabilities[row] * loopProbability; + } + } + + STORM_LOG_TRACE((hasSelfLoop ? "State has self-loop." : "State does not have a self-loop.")); + + // Now connect the predecessors of the state being eliminated with its successors. + row_type& currentStatePredecessors = backwardTransitions.getRow(state); + std::size_t predecessorForwardTransitionCount = 0; + + // In case we have a constrained elimination, we need to keep track of the new predecessors. + row_type newCurrentStatePredecessors; + + // Now go through the predecessors and eliminate the ones (satisfying the constraint if given). + for (auto const& predecessorEntry : currentStatePredecessors) { + uint_fast64_t predecessor = predecessorEntry.getColumn(); + + // Skip the state itself as one of its predecessors. + if (predecessor == row) { + assert(hasSelfLoop); + continue; + } + + // Skip the state if the elimination is constrained, but the predecessor is not in the constraint. + if (constrained && !predecessorConstraint.get(predecessor)) { + newCurrentStatePredecessors.emplace_back(predecessor, storm::utility::one()); + STORM_LOG_TRACE("Not eliminating predecessor " << predecessor << ", because it does not fit the filter."); + continue; + } + STORM_LOG_TRACE("Eliminating predecessor " << predecessor << "."); + + // First, find the probability with which the predecessor can move to the current state, because + // the other probabilities need to be scaled with this factor. + row_type& predecessorForwardTransitions = matrix.getRow(predecessor); + predecessorForwardTransitionCount += predecessorForwardTransitions.size(); + typename row_type::iterator multiplyElement = std::find_if(predecessorForwardTransitions.begin(), predecessorForwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == state; }); + + // Make sure we have found the probability and set it to zero. + STORM_LOG_THROW(multiplyElement != predecessorForwardTransitions.end(), storm::exceptions::InvalidStateException, "No probability for successor found."); + ValueType multiplyFactor = multiplyElement->getValue(); + multiplyElement->setValue(storm::utility::zero()); + + // At this point, we need to update the (forward) transitions of the predecessor. + typename row_type::iterator first1 = predecessorForwardTransitions.begin(); + typename row_type::iterator last1 = predecessorForwardTransitions.end(); + typename row_type::iterator first2 = currentStateSuccessors.begin(); + typename row_type::iterator last2 = currentStateSuccessors.end(); + + row_type newSuccessors; + newSuccessors.reserve((last1 - first1) + (last2 - first2)); + std::insert_iterator result(newSuccessors, newSuccessors.end()); + + // Now we merge the two successor lists. (Code taken from std::set_union and modified to suit our needs). + for (; first1 != last1; ++result) { + // Skip the transitions to the state that is currently being eliminated. + if (first1->getColumn() == state || (first2 != last2 && first2->getColumn() == state)) { + if (first1->getColumn() == state) { + ++first1; + } + if (first2 != last2 && first2->getColumn() == state) { + ++first2; + } + continue; + } + + if (first2 == last2) { + std::copy_if(first1, last1, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; } ); + break; + } + if (first2->getColumn() < first1->getColumn()) { + *result = storm::utility::simplify(std::move(*first2 * multiplyFactor)); + ++first2; + } else if (first1->getColumn() < first2->getColumn()) { + *result = *first1; + ++first1; + } else { + *result = storm::storage::MatrixEntry(first1->getColumn(), storm::utility::simplify(first1->getValue() + storm::utility::simplify(multiplyFactor * first2->getValue()))); + ++first1; + ++first2; + } + } + for (; first2 != last2; ++first2) { + if (first2->getColumn() != state) { + *result = storm::utility::simplify(std::move(*first2 * multiplyFactor)); + } + } + + // Now move the new transitions in place. + predecessorForwardTransitions = std::move(newSuccessors); + + if (!stateRewards) { + // Add the probabilities to go to a target state in just one step if we have to compute probabilities. + oneStepProbabilities[predecessor] += storm::utility::simplify(multiplyFactor * oneStepProbabilities[row]); + STORM_LOG_TRACE("Fixed new next-state probabilities of predecessor states."); + } else { + // If we are computing rewards, we basically scale the state reward of the state to eliminate and + // add the result to the state reward of the predecessor. + if (hasSelfLoop) { + stateRewards.get()[predecessor] += storm::utility::simplify(multiplyFactor * loopProbability * stateRewards.get()[row]); + } else { + stateRewards.get()[predecessor] += storm::utility::simplify(multiplyFactor * stateRewards.get()[row]); + } + } + } + + // Finally, we need to add the predecessor to the set of predecessors of every successor. + for (auto const& successorEntry : currentStateSuccessors) { + row_type& successorBackwardTransitions = backwardTransitions.getRow(successorEntry.getColumn()); + + // Delete the current state as a predecessor of the successor state only if we are going to remove the + // current state's forward transitions. + if (removeForwardTransitions) { + typename row_type::iterator elimIt = std::find_if(successorBackwardTransitions.begin(), successorBackwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == row; }); + STORM_LOG_ASSERT(elimIt != successorBackwardTransitions.end(), "Expected a proper backward transition, but found none."); + successorBackwardTransitions.erase(elimIt); + } + + typename row_type::iterator first1 = successorBackwardTransitions.begin(); + typename row_type::iterator last1 = successorBackwardTransitions.end(); + typename row_type::iterator first2 = currentStatePredecessors.begin(); + typename row_type::iterator last2 = currentStatePredecessors.end(); + + row_type newPredecessors; + newPredecessors.reserve((last1 - first1) + (last2 - first2)); + std::insert_iterator result(newPredecessors, newPredecessors.end()); + + if (!constrained) { + for (; first1 != last1; ++result) { + if (first2 == last2) { + std::copy(first1, last1, result); + break; + } + if (first2->getColumn() < first1->getColumn()) { + if (first2->getColumn() != row) { + *result = *first2; + } + ++first2; + } else { + *result = *first1; + if (first1->getColumn() == first2->getColumn()) { + ++first2; + } + ++first1; + } + } + std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != row; }); + } else { + // If the elimination is constrained, we need to be more selective when we set the new predecessors + // of the successor state. + for (; first1 != last1; ++result) { + if (first2 == last2) { + std::copy(first1, last1, result); + break; + } + if (first2->getColumn() < first1->getColumn()) { + if (first2->getColumn() != row) { + *result = *first2; + } + ++first2; + } else { + *result = *first1; + if (first1->getColumn() == first2->getColumn()) { + ++first2; + } + ++first1; + } + } + std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != row && (!constrained || predecessorConstraint.get(a.getColumn())); }); + } + + // Now move the new predecessors in place. + successorBackwardTransitions = std::move(newPredecessors); + } + STORM_LOG_TRACE("Fixed predecessor lists of successor states."); + + if (removeForwardTransitions) { + // Clear the eliminated row to reduce memory consumption. + currentStateSuccessors.clear(); + currentStateSuccessors.shrink_to_fit(); + } + if (!constrained) { + currentStatePredecessors.clear(); + currentStatePredecessors.shrink_to_fit(); + } else { + currentStatePredecessors = std::move(newCurrentStatePredecessors); + } + } + + template class FlexibleSparseMatrix; + template class FlexibleSparseMatrix; +#ifdef STORM_HAVE_CARL + template class FlexibleSparseMatrix; +#endif + + + } +} \ No newline at end of file diff --git a/src/storage/FlexibleSparseMatrix.h b/src/storage/FlexibleSparseMatrix.h new file mode 100644 index 000000000..d29b46995 --- /dev/null +++ b/src/storage/FlexibleSparseMatrix.h @@ -0,0 +1,51 @@ +#ifndef STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ +#define STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ + +#include + +#include "src/storage/SparseMatrix.h" +#include "src/storage/sparse/StateType.h" +#include "src/storage/BitVector.h" + +namespace storm { + namespace storage { + + template + class FlexibleSparseMatrix { + public: + typedef uint_fast64_t index_type; + typedef std::vector> row_type; + typedef typename row_type::iterator iterator; + typedef typename row_type::const_iterator const_iterator; + + FlexibleSparseMatrix() = default; + FlexibleSparseMatrix(index_type rows); + FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false); + + void reserveInRow(index_type row, index_type numberOfElements); + + row_type& getRow(index_type); + row_type const& getRow(index_type) const; + + index_type getNumberOfRows() const; + + void print() const; + + /*! + * Checks whether the given state has a self-loop with an arbitrary probability in the given probability matrix. + * + * @param state The state for which to check whether it possesses a self-loop. + * @param matrix The matrix in which to look for the loop. + * @return True iff the given state has a self-loop with an arbitrary probability in the given probability matrix. + */ + bool hasSelfLoop(storm::storage::sparse::state_type state); + + static void eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, uint_fast64_t row, FlexibleSparseMatrix& backwardTransitions, boost::optional>& stateRewards, bool removeForwardTransitions = true, bool constrained = false, storm::storage::BitVector const& predecessorConstraint = storm::storage::BitVector()); + + private: + std::vector data; + }; + } +} +#endif /* STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ */ + diff --git a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp index a4026bfe1..ab0a4e03d 100644 --- a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp @@ -98,7 +98,7 @@ TEST(SparseMdpRegionModelCheckerTest, coin_Prob) { EXPECT_NEAR(0.26787251126, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.41879628383, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.01535089684, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); - EXPECT_NEAR(0.24952612245, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.24952471590, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.01711494956, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); //test approximative method