diff --git a/src/modelchecker/reachability/SparseSccModelChecker.cpp b/src/modelchecker/reachability/SparseSccModelChecker.cpp index d9226ad66..78c54a4a0 100644 --- a/src/modelchecker/reachability/SparseSccModelChecker.cpp +++ b/src/modelchecker/reachability/SparseSccModelChecker.cpp @@ -4,6 +4,8 @@ #include "src/storage/StronglyConnectedComponentDecomposition.h" #include "src/utility/graph.h" +#include "src/utility/vector.h" +#include "src/exceptions/InvalidStateException.h" #include "src/exceptions/ExceptionMacros.h" namespace storm { @@ -30,41 +32,35 @@ namespace storm { // 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(dtmc.getTransitionMatrix(), dtmc.getInitialStates(), maybeStates, statesWithProbability1); - // Subtract from the maybe states the set of states that is not reachable. + // 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; // Otherwise, we build the submatrix that only has the transitions of the maybe states. storm::storage::SparseMatrix submatrix = dtmc.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates); // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. - FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(dtmc.getTransitionMatrix()); + FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(submatrix); // Create a vector for the probabilities to go to a state with probability 1 in one step. std::vector oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); - for (uint_fast64_t i = 0; i < oneStepProbabilities.size(); ++i) { - std::cout << i << " -> " << oneStepProbabilities[i] << std::endl; - } // Then, we recursively treat all SCCs. - treatScc(dtmc, flexibleMatrix, oneStepProbabilities, dtmc.getInitialStates() % maybeStates, storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true), submatrix, submatrix.transpose(), false, 0); + FlexibleSparseMatrix backwardTransitions = getFlexibleSparseMatrix(submatrix.transpose(), true); + treatScc(dtmc, flexibleMatrix, oneStepProbabilities, dtmc.getInitialStates() % maybeStates, storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true), submatrix, backwardTransitions, false, 0); // Now, we return the value for the only initial state. return oneStepProbabilities[initialStateIndex]; } template - void SparseSccModelChecker::treatScc(storm::models::Dtmc const& dtmc, FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, storm::storage::SparseMatrix const& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level) { - std::cout << "treating SCC " << scc << " at level " << level << " (entry: " << entryStates << ")" << std::endl; - if (level <= 2) { - std::cout << "1" << std::endl; + void SparseSccModelChecker::treatScc(storm::models::Dtmc const& dtmc, 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) { + if (level <= 3) { // Here, we further decompose the SCC into sub-SCCs. storm::storage::StronglyConnectedComponentDecomposition decomposition(forwardTransitions, scc & ~entryStates, true, false); - std::cout << decomposition << std::endl; // To eliminate the remaining one-state SCCs, we need to keep track of them. storm::storage::BitVector remainingStates(scc); - std::cout << "1" << std::endl; // And then recursively treat all sub-SCCs. for (auto const& newScc : decomposition) { // If the SCC consists of just one state, we do not explore it recursively, but rather eliminate @@ -90,8 +86,6 @@ namespace storm { // Recursively descend in SCC-hierarchy. treatScc(dtmc, matrix, oneStepProbabilities, entryStates, newSccAsBitVector, forwardTransitions, backwardTransitions, true, level + 1); } - - std::cout << "1" << std::endl; // If we are not supposed to eliminate the entry states, we need to take them out of the set of // remaining states. @@ -122,13 +116,16 @@ namespace storm { } template - void SparseSccModelChecker::eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, storm::storage::SparseMatrix const& backwardTransitions) { + void SparseSccModelChecker::eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix& backwardTransitions) { ValueType loopProbability = storm::utility::constantZero(); - + // Start by finding loop probability. - for (auto const& entry : matrix.getRow(state)) { - if (entry.getColumn() == state) { - loopProbability = entry.getValue(); + typename FlexibleSparseMatrix::row_type& currentStateSuccessors = matrix.getRow(state); + for (auto const& entry : currentStateSuccessors) { + if (entry.getColumn() >= state) { + if (entry.getColumn() == state) { + loopProbability = entry.getValue(); + } break; } } @@ -140,49 +137,128 @@ namespace storm { } oneStepProbabilities[state] *= loopProbability; - // Now connect the predecessors of the state to eliminate with its successors. - std::size_t newEntries = matrix.getRow(state).size(); - for (auto const& predecessorEntry : backwardTransitions.getRow(state)) { - // First, add all entries of the successor to the list of outgoing transitions of the predecessor. - typename FlexibleSparseMatrix::row_type row = matrix.getRow(predecessorEntry.getColumn()); - typename FlexibleSparseMatrix::row_type::iterator multiplyElement = std::find_if(row.begin(), row.end(), [=](storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() == state; }); - - ValueType multiplyFactor = storm::utility::constantOne(); - if (multiplyElement != row.end()) { - // Remove the transition to the state that is to be eliminated. - multiplyElement->setValue(0); - multiplyFactor = multiplyElement->getValue(); + // Now connect the predecessors of the state being eliminated with its successors. + typename FlexibleSparseMatrix::row_type& currentStatePredecessors = backwardTransitions.getRow(state); + for (auto const& predecessorEntry : currentStatePredecessors) { + uint_fast64_t predecessor = predecessorEntry.getColumn(); + + // Skip the state itself as one of its predecessors. + if (predecessor == state) { + continue; } - // Now scale all the entries in the current row and insert them in the transitions of the predecessor. - row.reserve(row.size() + newEntries); - std::for_each(matrix.getRow(state).begin(), matrix.getRow(state).end(), [&] (storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { row.emplace_back(a.getColumn(), multiplyFactor * a.getValue()); }); - - // Then sort the vector according to their column indices. - std::sort(row.begin(), row.end(), [](storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a, storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& b){ return a.getColumn() < b.getColumn(); }); - - // Now we can eliminate entries with the same column by simple addition. - typename FlexibleSparseMatrix::row_type::iterator rowIt = row.begin(); - typename FlexibleSparseMatrix::row_type::iterator it = row.begin(); - typename FlexibleSparseMatrix::row_type::iterator rowIte = row.end(); - for (++it; it != rowIte; ++it) { - if (it->getValue() == storm::utility::constantZero()) { + // 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); + typename FlexibleSparseMatrix::row_type::iterator multiplyElement = std::find_if(predecessorForwardTransitions.begin(), predecessorForwardTransitions.end(), [&](storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() == state; }); + + // Make sure we have found the probability and set it to zero. + LOG_THROW(multiplyElement != predecessorForwardTransitions.end(), storm::exceptions::InvalidStateException, "No probability for successor found."); + ValueType multiplyFactor = multiplyElement->getValue(); + multiplyElement->setValue(0); + + // 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::row_type> 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 (it->getColumn() == rowIt->getColumn()) { - rowIt->setValue(rowIt->getValue() + it->getValue()); + if (first2 == last2) { + std::copy_if(first1, last1, result, [&] (storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() != state; } ); + break; + } + if (first2->getColumn() < first1->getColumn()) { + *result = *first2 * multiplyFactor; + ++first2; + } else if (first1->getColumn() < first2->getColumn()) { + *result = *first1; + ++first1; } else { - ++rowIt; - *rowIt = *it; + *result = storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type>(first1->getColumn(), first1->getValue() + multiplyFactor * first2->getValue()); + ++first1; + ++first2; } } + for (; first2 != last2; ++first2) { + if (first2->getColumn() != state) { + *result = *first2 * multiplyFactor; + } + } + + // Now move the new transitions in place. + predecessorForwardTransitions = std::move(newSuccessors); - // Finally, add the probabilities to go to a target state in just one step. - std::cout << "prior: " << oneStepProbabilities[predecessorEntry.getColumn()] << std::endl; - oneStepProbabilities[predecessorEntry.getColumn()] += oneStepProbabilities[state]; - std::cout << "updating one step prob of " << predecessorEntry.getColumn() << " to " << oneStepProbabilities[predecessorEntry.getColumn()] << std::endl; + // Add the probabilities to go to a target state in just one step. + oneStepProbabilities[predecessor] += multiplyFactor * oneStepProbabilities[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. + typename FlexibleSparseMatrix::row_type::const_iterator elimIt = std::find_if(successorBackwardTransitions.begin(), successorBackwardTransitions.end(), [&](storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() == state; }); + if (elimIt != successorBackwardTransitions.end()) { + 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::row_type> result(newPredecessors, newPredecessors.end()); + + + for (; first1 != last1; ++result) { + if (first2 == last2) { + std::copy_if(first1, last1, result, [&] (storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() != state; }); + break; + } + if (first2->getColumn() < first1->getColumn()) { + if (first2->getColumn() != state) { + *result = *first2; + } + ++first2; + } else { + if (first1->getColumn() != state) { + *result = *first1; + } + if (first1->getColumn() == first2->getColumn()) { + ++first2; + } + ++first1; + } + } + std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type> const& a) { return a.getColumn() != state; }); + + // Now move the new predecessors in place. + successorBackwardTransitions = std::move(newPredecessors); + } + + + // Clear the eliminated row to reduce memory consumption. + currentStateSuccessors.clear(); + currentStateSuccessors.shrink_to_fit(); } template @@ -201,7 +277,7 @@ namespace storm { } template - FlexibleSparseMatrix SparseSccModelChecker::getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix) { + FlexibleSparseMatrix SparseSccModelChecker::getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne) { FlexibleSparseMatrix flexibleMatrix(matrix.getRowCount()); for (typename FlexibleSparseMatrix::index_type rowIndex = 0; rowIndex < matrix.getRowCount(); ++rowIndex) { @@ -209,7 +285,11 @@ namespace storm { flexibleMatrix.reserveInRow(rowIndex, row.getNumberOfEntries()); for (auto const& element : row) { - flexibleMatrix.getRow(rowIndex).emplace_back(element); + if (setAllValuesToOne) { + flexibleMatrix.getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::constantOne()); + } else { + flexibleMatrix.getRow(rowIndex).emplace_back(element); + } } } diff --git a/src/modelchecker/reachability/SparseSccModelChecker.h b/src/modelchecker/reachability/SparseSccModelChecker.h index ecb97e8b1..ae3ff102e 100644 --- a/src/modelchecker/reachability/SparseSccModelChecker.h +++ b/src/modelchecker/reachability/SparseSccModelChecker.h @@ -33,9 +33,9 @@ namespace storm { static ValueType computeReachabilityProbability(storm::models::Dtmc const& dtmc, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates); private: - static void treatScc(storm::models::Dtmc const& dtmc, FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, storm::storage::SparseMatrix const& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level); - static FlexibleSparseMatrix getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix); - static void eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, storm::storage::SparseMatrix const& backwardTransitions); + static void treatScc(storm::models::Dtmc const& dtmc, 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); + static FlexibleSparseMatrix getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false); + static void eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix& backwardTransitions); }; } diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 010cc36bf..e52bd87ef 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -53,6 +53,17 @@ namespace storm { return this->entry; } + template + MatrixEntry MatrixEntry::operator*(value_type factor) const { + return MatrixEntry(this->getColumn(), this->getValue() * factor); + } + + template + std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry) { + out << "(" << entry.getColumn() << ", " << entry.getValue() << ")"; + return out; + } + template SparseMatrixBuilder::SparseMatrixBuilder(index_type rows, index_type columns, index_type entries, bool forceDimensions, bool hasCustomRowGrouping, index_type rowGroups) : initialRowCountSet(rows != 0), initialRowCount(rows), initialColumnCountSet(columns != 0), initialColumnCount(columns), initialEntryCountSet(entries != 0), initialEntryCount(entries), forceInitialDimensions(forceDimensions), hasCustomRowGrouping(hasCustomRowGrouping), initialRowGroupCountSet(rowGroups != 0), initialRowGroupCount(rowGroups), rowGroupIndices(), columnsAndValues(), rowIndications(), currentEntryCount(0), lastRow(0), lastColumn(0), highestColumn(0), currentRowGroup(0) { // Prepare the internal storage. @@ -948,6 +959,7 @@ namespace storm { // Explicitly instantiate the entry, builder and the matrix. template class MatrixEntry::index_type, double>; + template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index b419ca841..856261c4f 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -95,6 +95,15 @@ namespace storm { */ std::pair const& getColumnValuePair() const; + /*! + * Multiplies the entry with the given factor and returns the result. + * + * @param factor The factor with which to multiply the entry. + */ + MatrixEntry operator*(value_type factor) const; + + template + friend std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); private: // The actual matrix entry. std::pair entry; diff --git a/src/storm.cpp b/src/storm.cpp index ff2e4da47..75ea53fb8 100644 --- a/src/storm.cpp +++ b/src/storm.cpp @@ -102,7 +102,7 @@ int main(const int argc, const char* argv[]) { initializeLogger(); setUp(); - try { +// try { LOG4CPLUS_INFO(logger, "StoRM was invoked."); // Parse options. @@ -229,10 +229,11 @@ int main(const int argc, const char* argv[]) { std::cout << "Parsing and translating the Symbolic Input took " << std::chrono::duration_cast(programTranslationEnd - programTranslationStart).count() << " milliseconds." << std::endl; storm::modelchecker::reachability::SparseSccModelChecker modelChecker; - storm::storage::BitVector trueStates(model->getNumberOfStates(), true); - storm::storage::BitVector oneStates = model->getLabeledStates("one"); - double value = modelChecker.computeReachabilityProbability(*model->as>(), trueStates, oneStates); + storm::storage::BitVector targetStates = model->getLabeledStates("observe0Greater1"); +// storm::storage::BitVector targetStates = model->getLabeledStates("one"); +// storm::storage::BitVector targetStates = model->getLabeledStates("elected"); + double value = modelChecker.computeReachabilityProbability(*model->as>(), trueStates, targetStates); std::cout << "computed value " << value << std::endl; if (s->isSet("mincmd")) { @@ -319,9 +320,9 @@ int main(const int argc, const char* argv[]) { printUsage(); LOG4CPLUS_INFO(logger, "StoRM terminating."); return 0; - } catch (std::exception& e) { - LOG4CPLUS_FATAL(logger, "An exception was thrown. Terminating."); - LOG4CPLUS_FATAL(logger, "\t" << e.what()); - } +// } catch (std::exception& e) { +// LOG4CPLUS_FATAL(logger, "An exception was thrown. Terminating."); +// LOG4CPLUS_FATAL(logger, "\t" << e.what()); +// } return 1; }