first version of conditional probabilities for (non-parametric) DTMCs a la Baier

Former-commit-id: b57dfab024

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@@ -33,7 +33,19 @@ namespace storm {
         
         template<typename SparseDtmcModelType>
         bool SparseDtmcPrctlModelChecker<SparseDtmcModelType>::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<typename SparseDtmcModelType>
@@ -95,7 +107,21 @@ namespace storm {
             std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeLongRunAverageProbabilities(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), nullptr, qualitative, *linearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
-                
+        
+        template<typename SparseDtmcModelType>
+        std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative , boost::optional<OptimizationDirection> 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<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula().asEventuallyFormula().getSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula().asEventuallyFormula().getSubformula());
+            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
+            ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
+
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeConditionalProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
+        }
+        
         template class SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>>;
     }
 }

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@@ -23,6 +23,7 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
+            virtual std::unique_ptr<CheckResult> computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
             virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
             virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;
             virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<OptimizationDirection> const& optimalityType = boost::optional<OptimizationDirection>()) override;

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<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return SparseCtmcCslHelper<ValueType>::computeLongRunAverageProbabilities(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory);
             }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                std::vector<ValueType> 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<ValueType>()) {
+                        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<uint_fast64_t> numberOfNormalStatesUpToState = statesWithProbabilityGreater0.getNumberOfSetBitsBeforeIndices();
+
+                // Now, we create the matrix of 'before' and 'normal' states.
+                std::vector<uint_fast64_t> numberOfBeforeStatesUpToState = beforeStates.getNumberOfSetBitsBeforeIndices();
+                storm::storage::SparseMatrixBuilder<ValueType> 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<ValueType>();
+                    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<ValueType>());
+
+                // Build the new transition matrix and the new targets.
+                storm::storage::SparseMatrix<ValueType> 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<ValueType> conditionalProbabilities = computeUntilProbabilities(newTransitionMatrix, newTransitionMatrix.transpose(), storm::storage::BitVector(newTransitionMatrix.getRowCount(), true), newTargetStates, qualitative, linearEquationSolverFactory);
+
+                // Unpack and return result.
+                std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::infinity<ValueType>());
+                storm::utility::vector::setVectorValues(result, beforeStates, conditionalProbabilities);
+                return result;
+            }
 
             template class SparseDtmcPrctlHelper<double>;
         }

src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h

@@ -12,6 +12,8 @@
 
 namespace storm {
     namespace modelchecker {
+        class CheckResult;
+        
         namespace helper {
             
             template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
@@ -33,6 +35,8 @@ namespace storm {
                 
                 static std::vector<ValueType> computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
+                static std::vector<ValueType> computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
             private:
                 static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             };

src/storage/BitVector.cpp

@@ -564,6 +564,21 @@ namespace storm {
 
             return result;
         }
+        
+        std::vector<uint_fast64_t> BitVector::getNumberOfSetBitsBeforeIndices() const {
+            std::vector<uint_fast64_t> 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<size_t> (bitCount);

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<uint_fast64_t> getNumberOfSetBitsBeforeIndices() const;
+                        
             /*!
              * Retrieves the number of bits this bit vector can store.
              *

src/storage/SparseMatrix.cpp

@@ -412,7 +412,6 @@ namespace storm {
         
         template<typename ValueType>
         void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
-            //SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
             this->nonzeroEntryCount = 0;
             for (auto const& element : *this) {
                 if (element.getValue() != storm::utility::zero<ValueType>()) {
@@ -426,6 +425,18 @@ namespace storm {
             this->nonzeroEntryCount += difference;
         }
         
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::updateDimensions() const {
+            this->nonzeroEntryCount = 0;
+            this->columnCount = 0;
+            for (auto const& element : *this) {
+                if (element.getValue() != storm::utility::zero<ValueType>()) {
+                    ++this->nonzeroEntryCount;
+                    this->columnCount = std::max(element.getColumn() + 1, this->columnCount);
+                }
+            }
+        }
+        
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::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<index_type> 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<index_type> columnBitsSetBeforeIndex = columnConstraint.getNumberOfSetBitsBeforeIndices();
             std::vector<index_type>* rowBitsSetBeforeIndex;
             if (&rowGroupConstraint == &columnConstraint) {
                 rowBitsSetBeforeIndex = &columnBitsSetBeforeIndex;
             } else {
-                rowBitsSetBeforeIndex = new std::vector<index_type>(rowCount);
-                
-                lastIndex = 0;
-                currentNumberOfSetBits = 0;
-                for (auto index : rowGroupConstraint) {
-                    while (lastIndex <= index) {
-                        rowBitsSetBeforeIndex->push_back(currentNumberOfSetBits);
-                        ++lastIndex;
-                    }
-                    ++currentNumberOfSetBits;
-                }
+                rowBitsSetBeforeIndex = new std::vector<index_type>(rowGroupConstraint.getNumberOfSetBitsBeforeIndices());
             }
             
             // Then, we need to determine the number of entries and the number of rows of the submatrix.

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;