From 33757633c83f942961bdb05ac1de19e1f8f37130 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sat, 26 Dec 2015 12:45:26 +0100 Subject: [PATCH] first version of conditional probabilities for (non-parametric) DTMCs a la Baier Former-commit-id: b57dfab0240611429d9c5c3f35acc272b1e51654 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 30 ++++++- .../prctl/SparseDtmcPrctlModelChecker.h | 1 + .../prctl/helper/SparseDtmcPrctlHelper.cpp | 88 ++++++++++++++++++- .../prctl/helper/SparseDtmcPrctlHelper.h | 4 + src/storage/BitVector.cpp | 15 ++++ src/storage/BitVector.h | 7 ++ src/storage/SparseMatrix.cpp | 39 +++----- src/storage/SparseMatrix.h | 7 +- 8 files changed, 162 insertions(+), 29 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index b4fa14522..9a7ae0899 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -33,7 +33,19 @@ namespace storm { template bool SparseDtmcPrctlModelChecker::canHandle(storm::logic::Formula const& formula) const { - return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula(); + if (formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula()) { + return true; + } + if (formula.isProbabilityOperatorFormula()) { + return this->canHandle(formula.asProbabilityOperatorFormula().getSubformula()); + } + if (formula.isConditionalPathFormula()) { + storm::logic::ConditionalPathFormula const& conditionalPathFormula = formula.asConditionalPathFormula(); + if (conditionalPathFormula.getLeftSubformula().isEventuallyFormula() && conditionalPathFormula.getRightSubformula().isEventuallyFormula()) { + return this->canHandle(conditionalPathFormula.getLeftSubformula()) && this->canHandle(conditionalPathFormula.getRightSubformula()); + } + } + return false; } template @@ -95,7 +107,21 @@ namespace storm { std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), nullptr, qualitative, *linearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } - + + template + std::unique_ptr SparseDtmcPrctlModelChecker::computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative , boost::optional const& optimalityType) { + STORM_LOG_THROW(pathFormula.getLeftSubformula().isEventuallyFormula(), storm::exceptions::InvalidPropertyException, "Illegal conditional probability formula."); + STORM_LOG_THROW(pathFormula.getRightSubformula().isEventuallyFormula(), storm::exceptions::InvalidPropertyException, "Illegal conditional probability formula."); + + std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula().asEventuallyFormula().getSubformula()); + std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula().asEventuallyFormula().getSubformula()); + ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); + ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); + + std::vector numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeConditionalProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + } + template class SparseDtmcPrctlModelChecker>; } } \ No newline at end of file diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index 14d888aa9..956ecceaa 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -23,6 +23,7 @@ namespace storm { virtual std::unique_ptr computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName = boost::optional(), bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName = boost::optional(), bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName = boost::optional(), bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 9c8f85385..bef7ede14 100644 --- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -6,9 +6,10 @@ #include "src/utility/vector.h" #include "src/utility/graph.h" - #include "src/solver/LinearEquationSolver.h" +#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" + #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" @@ -217,6 +218,91 @@ namespace storm { std::vector SparseDtmcPrctlHelper::computeLongRunAverageProbabilities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { return SparseCtmcCslHelper::computeLongRunAverageProbabilities(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory); } + + template + std::vector SparseDtmcPrctlHelper::computeConditionalProbabilities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + std::vector probabilitiesToReachConditionStates = computeUntilProbabilities(transitionMatrix, backwardTransitions, storm::storage::BitVector(transitionMatrix.getRowCount(), true), conditionStates, false, linearEquationSolverFactory); + + // Start by computing all 'before' states, i.e. the states for which the conditional probability is defined. + storm::storage::BitVector beforeStates(targetStates.size(), true); + uint_fast64_t state = 0; + uint_fast64_t beforeStateIndex = 0; + for (auto const& value : probabilitiesToReachConditionStates) { + if (value == storm::utility::zero()) { + beforeStates.set(state, false); + } else { + probabilitiesToReachConditionStates[beforeStateIndex] = value; + ++beforeStateIndex; + } + ++state; + } + probabilitiesToReachConditionStates.resize(beforeStateIndex); + + uint_fast64_t normalStatesOffset = beforeStates.getNumberOfSetBits(); + storm::storage::BitVector allStates(targetStates.size(), true); + storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(backwardTransitions, allStates, targetStates); + statesWithProbabilityGreater0 &= storm::utility::graph::getReachableStates(transitionMatrix, conditionStates, allStates, targetStates); + std::vector numberOfNormalStatesUpToState = statesWithProbabilityGreater0.getNumberOfSetBitsBeforeIndices(); + + // Now, we create the matrix of 'before' and 'normal' states. + std::vector numberOfBeforeStatesUpToState = beforeStates.getNumberOfSetBitsBeforeIndices(); + storm::storage::SparseMatrixBuilder builder; + uint_fast64_t currentRow = 0; + + // Start by creating the transitions of the 'before' states. + for (auto beforeState : beforeStates) { + if (conditionStates.get(beforeState)) { + // For condition states, we move to the 'normal' states. + for (auto const& successorEntry : transitionMatrix.getRow(beforeState)) { + if (statesWithProbabilityGreater0.get(successorEntry.getColumn())) { + builder.addNextValue(currentRow, normalStatesOffset + numberOfNormalStatesUpToState[successorEntry.getColumn()], successorEntry.getValue()); + } + } + } else { + // For non-condition states, we scale the probabilities going to other before states. + for (auto const& successorEntry : transitionMatrix.getRow(beforeState)) { + if (beforeStates.get(successorEntry.getColumn())) { + builder.addNextValue(currentRow, numberOfBeforeStatesUpToState[successorEntry.getColumn()], successorEntry.getValue() * probabilitiesToReachConditionStates[numberOfBeforeStatesUpToState[successorEntry.getColumn()]] / probabilitiesToReachConditionStates[currentRow]); + } + } + } + ++currentRow; + } + + // Then, create the transitions of the 'normal' states. + uint_fast64_t deadlockState = normalStatesOffset + statesWithProbabilityGreater0.getNumberOfSetBits(); + for (auto state : statesWithProbabilityGreater0) { + ValueType zeroProbability = storm::utility::zero(); + for (auto const& successorEntry : transitionMatrix.getRow(state)) { + if (statesWithProbabilityGreater0.get(successorEntry.getColumn())) { + builder.addNextValue(currentRow, normalStatesOffset + numberOfNormalStatesUpToState[successorEntry.getColumn()], successorEntry.getValue()); + } else { + zeroProbability += successorEntry.getValue(); + } + } + if (!storm::utility::isZero(zeroProbability)) { + builder.addNextValue(currentRow, deadlockState, zeroProbability); + } + ++currentRow; + } + builder.addNextValue(deadlockState, deadlockState, storm::utility::one()); + + // Build the new transition matrix and the new targets. + storm::storage::SparseMatrix newTransitionMatrix = builder.build(); + storm::storage::BitVector newTargetStates = targetStates % beforeStates; + newTargetStates.resize(newTransitionMatrix.getRowCount()); + for (auto state : targetStates % statesWithProbabilityGreater0) { + newTargetStates.set(normalStatesOffset + state, true); + } + + // Finally, compute the conditional probabiltities by a reachability query. + std::vector conditionalProbabilities = computeUntilProbabilities(newTransitionMatrix, newTransitionMatrix.transpose(), storm::storage::BitVector(newTransitionMatrix.getRowCount(), true), newTargetStates, qualitative, linearEquationSolverFactory); + + // Unpack and return result. + std::vector result(transitionMatrix.getRowCount(), storm::utility::infinity()); + storm::utility::vector::setVectorValues(result, beforeStates, conditionalProbabilities); + return result; + } template class SparseDtmcPrctlHelper; } diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h index c28fed3b8..ff83dac63 100644 --- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h +++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h @@ -12,6 +12,8 @@ namespace storm { namespace modelchecker { + class CheckResult; + namespace helper { template > @@ -33,6 +35,8 @@ namespace storm { static std::vector computeLongRunAverageProbabilities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + static std::vector computeConditionalProbabilities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + private: static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); }; diff --git a/src/storage/BitVector.cpp b/src/storage/BitVector.cpp index 0dcc3cbec..cd5f46009 100644 --- a/src/storage/BitVector.cpp +++ b/src/storage/BitVector.cpp @@ -564,6 +564,21 @@ namespace storm { return result; } + + std::vector BitVector::getNumberOfSetBitsBeforeIndices() const { + std::vector bitsSetBeforeIndices; + bitsSetBeforeIndices.reserve(this->size()); + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + for (auto index : *this) { + while (lastIndex <= index) { + bitsSetBeforeIndices.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + return bitsSetBeforeIndices; + } size_t BitVector::size() const { return static_cast (bitCount); diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index 8aed9e387..f9c088ceb 100644 --- a/src/storage/BitVector.h +++ b/src/storage/BitVector.h @@ -412,6 +412,13 @@ namespace storm { */ uint_fast64_t getNumberOfSetBitsBeforeIndex(uint_fast64_t index) const; + /*! + * Retrieves a vector that holds at position i the number of bits set before index i. + * + * @return The resulting vector of 'offsets'. + */ + std::vector getNumberOfSetBitsBeforeIndices() const; + /*! * Retrieves the number of bits this bit vector can store. * diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 6850eee09..5130ffc5f 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -412,7 +412,6 @@ namespace storm { template void SparseMatrix::updateNonzeroEntryCount() const { - //SparseMatrix* self = const_cast*>(this); this->nonzeroEntryCount = 0; for (auto const& element : *this) { if (element.getValue() != storm::utility::zero()) { @@ -426,6 +425,18 @@ namespace storm { this->nonzeroEntryCount += difference; } + template + void SparseMatrix::updateDimensions() const { + this->nonzeroEntryCount = 0; + this->columnCount = 0; + for (auto const& element : *this) { + if (element.getValue() != storm::utility::zero()) { + ++this->nonzeroEntryCount; + this->columnCount = std::max(element.getColumn() + 1, this->columnCount); + } + } + } + template typename SparseMatrix::index_type SparseMatrix::getRowGroupCount() const { return rowGroupIndices.size() - 1; @@ -554,34 +565,12 @@ namespace storm { // Start by creating a temporary vector that stores for each index whose bit is set to true the number of // bits that were set before that particular index. - std::vector columnBitsSetBeforeIndex; - columnBitsSetBeforeIndex.reserve(columnCount); - - // Compute the information to fill this vector. - index_type lastIndex = 0; - index_type currentNumberOfSetBits = 0; - for (auto index : columnConstraint) { - while (lastIndex <= index) { - columnBitsSetBeforeIndex.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } + std::vector columnBitsSetBeforeIndex = columnConstraint.getNumberOfSetBitsBeforeIndices(); std::vector* rowBitsSetBeforeIndex; if (&rowGroupConstraint == &columnConstraint) { rowBitsSetBeforeIndex = &columnBitsSetBeforeIndex; } else { - rowBitsSetBeforeIndex = new std::vector(rowCount); - - lastIndex = 0; - currentNumberOfSetBits = 0; - for (auto index : rowGroupConstraint) { - while (lastIndex <= index) { - rowBitsSetBeforeIndex->push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } + rowBitsSetBeforeIndex = new std::vector(rowGroupConstraint.getNumberOfSetBitsBeforeIndices()); } // Then, we need to determine the number of entries and the number of rows of the submatrix. diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index c31e0b344..b10172f65 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -514,6 +514,11 @@ namespace storm { * Recompute the nonzero entry count */ void updateNonzeroEntryCount() const; + + /*! + * Recomputes the number of columns and the number of non-zero entries. + */ + void updateDimensions() const; /*! * Change the nonzero entry count by the provided value. @@ -931,7 +936,7 @@ namespace storm { index_type rowCount; // The number of columns of the matrix. - index_type columnCount; + mutable index_type columnCount; // The number of entries in the matrix. index_type entryCount;