some more optimizations

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -960,6 +960,8 @@ namespace storm {
             template<typename ValueType>
             std::unique_ptr<CheckResult> SparseMdpPrctlHelper<ValueType>::computeConditionalProbabilities(OptimizationDirection dir, storm::storage::sparse::state_type initialState, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
                 
+                std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
+                
                 // For the max-case, we can simply take the given target states. For the min-case, however, we need to
                 // find the MECs of non-target states and make them the new target states.
                 storm::storage::BitVector fixedTargetStates;
@@ -990,7 +992,10 @@ namespace storm {
                 storm::storage::BitVector extendedConditionStates = storm::utility::graph::performProb1A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, allStates, conditionStates);
 
                 STORM_LOG_DEBUG("Computing probabilities to satisfy condition.");
+                std::chrono::high_resolution_clock::time_point conditionStart = std::chrono::high_resolution_clock::now();
                 std::vector<ValueType> conditionProbabilities = std::move(computeUntilProbabilities(OptimizationDirection::Maximize, transitionMatrix, backwardTransitions, allStates, extendedConditionStates, false, false, minMaxLinearEquationSolverFactory).values);
+                std::chrono::high_resolution_clock::time_point conditionEnd = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Computed probabilities to satisfy for condition in " << std::chrono::duration_cast<std::chrono::milliseconds>(conditionEnd - conditionStart).count() << "ms.");
                 
                 // If the conditional probability is undefined for the initial state, we return directly.
                 if (storm::utility::isZero(conditionProbabilities[initialState])) {
@@ -998,7 +1003,10 @@ namespace storm {
                 }
                 
                 STORM_LOG_DEBUG("Computing probabilities to reach target.");
+                std::chrono::high_resolution_clock::time_point targetStart = std::chrono::high_resolution_clock::now();
                 std::vector<ValueType> targetProbabilities = std::move(computeUntilProbabilities(OptimizationDirection::Maximize, transitionMatrix, backwardTransitions, allStates, fixedTargetStates, false, false, minMaxLinearEquationSolverFactory).values);
+                std::chrono::high_resolution_clock::time_point targetEnd = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Computed probabilities to reach target in " << std::chrono::duration_cast<std::chrono::milliseconds>(targetEnd - targetStart).count() << "ms.");
 
                 storm::storage::BitVector statesWithProbabilityGreater0E(transitionMatrix.getRowGroupCount(), true);
                 storm::storage::sparse::state_type state = 0;
@@ -1011,10 +1019,8 @@ namespace storm {
 
                 // Determine those states that need to be equipped with a restart mechanism.
                 STORM_LOG_DEBUG("Computing problematic states.");
-                storm::storage::BitVector pureResetStates = storm::utility::graph::performProb0A(backwardTransitions, allStates, fixedTargetStates);
-                
-                // FIXME: target | condition as target states here?
-                storm::storage::BitVector problematicStates = storm::utility::graph::performProb0E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, allStates, fixedTargetStates);
+                storm::storage::BitVector pureResetStates = storm::utility::graph::performProb0A(backwardTransitions, allStates, extendedConditionStates);
+                storm::storage::BitVector problematicStates = storm::utility::graph::performProb0E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, allStates, extendedConditionStates | fixedTargetStates);
 
                 // Otherwise, we build the transformed MDP.
                 storm::storage::BitVector relevantStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStatesBitVector, allStates, extendedConditionStates | fixedTargetStates | pureResetStates);
@@ -1073,6 +1079,9 @@ namespace storm {
                 builder.addNextValue(currentRow, numberOfStatesBeforeRelevantStates[initialState], storm::utility::one<ValueType>());
                 ++currentRow;
                 
+                std::chrono::high_resolution_clock::time_point end = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Computed transformed model in " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << "ms.");
+                
                 // Finally, build the matrix and dispatch the query as a reachability query.
                 STORM_LOG_DEBUG("Computing conditional probabilties.");
                 storm::storage::BitVector newGoalStates(newFailState + 1);
@@ -1080,7 +1089,10 @@ namespace storm {
                 storm::storage::SparseMatrix<ValueType> newTransitionMatrix = builder.build();
                 STORM_LOG_DEBUG("Transformed model has " << newTransitionMatrix.getRowGroupCount() << " states and " << newTransitionMatrix.getNonzeroEntryCount() << " transitions.");
                 storm::storage::SparseMatrix<ValueType> newBackwardTransitions = newTransitionMatrix.transpose(true);
+                std::chrono::high_resolution_clock::time_point conditionalStart = std::chrono::high_resolution_clock::now();
                 std::vector<ValueType> goalProbabilities = std::move(computeUntilProbabilities(OptimizationDirection::Maximize, newTransitionMatrix, newBackwardTransitions, storm::storage::BitVector(newFailState + 1, true), newGoalStates, false, false, minMaxLinearEquationSolverFactory).values);
+                std::chrono::high_resolution_clock::time_point conditionalEnd = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Computed conditional probabilities in transformed model in " << std::chrono::duration_cast<std::chrono::milliseconds>(conditionalEnd - conditionalStart).count() << "ms.");
                 
                 return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(initialState, dir == OptimizationDirection::Maximize ? goalProbabilities[numberOfStatesBeforeRelevantStates[initialState]] : storm::utility::one<ValueType>() - goalProbabilities[numberOfStatesBeforeRelevantStates[initialState]]));
             }

src/storm/models/sparse/StandardRewardModel.cpp

@@ -179,12 +179,12 @@ namespace storm {
                     STORM_LOG_THROW(transitionMatrix.getRowGroupCount() == this->getStateActionRewardVector().size(), storm::exceptions::InvalidOperationException, "The reduction to state rewards is only possible if the size of the action reward vector equals the number of states.");
                     if (weights) {
                         if (this->hasStateRewards()) {
-                            storm::utility::vector::applyPointwise<ValueType, MatrixValueType, ValueType>(this->getStateActionRewardVector(), *weights, this->getStateRewardVector(),
+                            storm::utility::vector::applyPointwiseTernary<ValueType, MatrixValueType, ValueType>(this->getStateActionRewardVector(), *weights, this->getStateRewardVector(),
                                                                    [] (ValueType const& sar, MatrixValueType const& w, ValueType const& sr) -> ValueType {
                                                                        return sr + w * sar; });
                         } else {
                             this->optionalStateRewardVector = std::move(this->optionalStateActionRewardVector);
-                            storm::utility::vector::applyPointwise<ValueType, MatrixValueType, ValueType>(this->optionalStateRewardVector.get(), *weights, this->optionalStateRewardVector.get(), [] (ValueType const& r, MatrixValueType const& w) { return w * r; } );
+                            storm::utility::vector::applyPointwise<ValueType, MatrixValueType, ValueType, std::multiplies<>>(this->optionalStateRewardVector.get(), *weights, this->optionalStateRewardVector.get());
                         }
                     } else {
                         if (this->hasStateRewards()) {
@@ -216,11 +216,11 @@ namespace storm {
                 std::vector<ValueType> result;
                 if (this->hasTransitionRewards()) {
                     result = transitionMatrix.getPointwiseProductRowSumVector(this->getTransitionRewardMatrix());
-                    storm::utility::vector::applyPointwise<MatrixValueType, ValueType, ValueType>(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * (resultElement + rewardElement); } );
+                    storm::utility::vector::applyPointwiseTernary<MatrixValueType, ValueType, ValueType>(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * (resultElement + rewardElement); } );
                 } else {
                     result = std::vector<ValueType>(transitionMatrix.getRowCount());
                     if (this->hasStateActionRewards()) {
-                        storm::utility::vector::applyPointwise<MatrixValueType, ValueType, ValueType>(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * rewardElement; } );
+                        storm::utility::vector::applyPointwise<MatrixValueType, ValueType, ValueType>(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement) { return weight * rewardElement; } );
                     }
                 }
                 if (this->hasStateRewards()) {

src/storm/utility/vector.h

@@ -307,8 +307,8 @@ namespace storm {
              * @param secondOperand The second operand.
              * @param target The target vector.
              */
-            template<class InValueType1, class InValueType2, class OutValueType>
-            void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, std::function<OutValueType (InValueType1 const&, InValueType2 const&, OutValueType const&)> const& function) {
+            template<class InValueType1, class InValueType2, class OutValueType, class Operation>
+            void applyPointwiseTernary(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, Operation f = Operation()) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
@@ -317,7 +317,7 @@ namespace storm {
                                       auto secondIt = secondOperand.begin() + range.begin();
                                       auto targetIt = target.begin() + range.begin();
                                       while (firstIt != firstIte) {
-                                          *targetIt = function(*firstIt, *secondIt, *targetIt);
+                                          *targetIt = f(*firstIt, *secondIt, *targetIt);
                                           ++targetIt;
                                           ++firstIt;
                                           ++secondIt;
@@ -329,7 +329,7 @@ namespace storm {
                 auto secondIt = secondOperand.begin();
                 auto targetIt = target.begin();
                 while (firstIt != firstIte) {
-                    *targetIt = function(*firstIt, *secondIt, *targetIt);
+                    *targetIt = f(*firstIt, *secondIt, *targetIt);
                     ++targetIt;
                     ++firstIt;
                     ++secondIt;
@@ -345,15 +345,15 @@ namespace storm {
              * @param secondOperand The second operand.
              * @param target The target vector.
              */
-            template<class InValueType1, class InValueType2, class OutValueType>
-            void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, std::function<OutValueType (InValueType1 const&, InValueType2 const&)> const& function) {
+            template<class InValueType1, class InValueType2, class OutValueType, class Operation>
+            void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, Operation f = Operation()) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
-                                      std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
+                                      std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), f);
                                   });
 #else
-                std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), function);
+                std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), f);
 #endif
             }
             
@@ -364,15 +364,15 @@ namespace storm {
              * @param target The target vector.
              * @param function The function to apply.
              */
-            template<class InValueType, class OutValueType>
-            void applyPointwise(std::vector<InValueType> const& operand, std::vector<OutValueType>& target, std::function<OutValueType (InValueType const&)> const& function) {
+            template<class InValueType, class OutValueType, class Operation>
+            void applyPointwise(std::vector<InValueType> const& operand, std::vector<OutValueType>& target, Operation f = Operation()) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
-                                      std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), function);
+                                      std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), f);
                                   });
 #else
-                std::transform(operand.begin(), operand.end(), target.begin(), function);
+                std::transform(operand.begin(), operand.end(), target.begin(), f);
 #endif
             }
             
@@ -385,7 +385,7 @@ namespace storm {
              */
             template<class InValueType1, class InValueType2, class OutValueType>
             void addVectors(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
-                applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::plus<>());
+                applyPointwise<InValueType1, InValueType2, OutValueType, std::plus<>>(firstOperand, secondOperand, target);
             }
             
             /*!
@@ -397,7 +397,7 @@ namespace storm {
              */
             template<class InValueType1, class InValueType2, class OutValueType>
             void subtractVectors(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
-                applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::minus<>());
+                applyPointwise<InValueType1, InValueType2, OutValueType, std::minus<>>(firstOperand, secondOperand, target);
             }
             
             /*!
@@ -409,7 +409,7 @@ namespace storm {
              */
             template<class InValueType1, class InValueType2, class OutValueType>
             void multiplyVectorsPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
-                applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::multiplies<>());
+                applyPointwise<InValueType1, InValueType2, OutValueType, std::multiplies<>>(firstOperand, secondOperand, target);
             }
             
             /*!
@@ -421,7 +421,7 @@ namespace storm {
              */
             template<class InValueType1, class InValueType2, class OutValueType>
             void divideVectorsPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
-                applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::divides<>());
+                applyPointwise<InValueType1, InValueType2, OutValueType, std::divides<>>(firstOperand, secondOperand, target);
             }
             
             /*!
@@ -603,8 +603,9 @@ namespace storm {
              * return true iff v1 is supposed to be taken instead of v2.
              * @param choices If non-null, this vector is used to store the choices made during the selection.
              */
-            template<class T>
-            void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::function<bool (T const&, T const&)> filter, std::vector<uint_fast64_t>* choices) {
+            template<class T, class Filter>
+            void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices) {
+                Filter f;
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
@@ -621,7 +622,7 @@ namespace storm {
                                       if (choices != nullptr) {
                                           choiceIt = choices->begin();
                                       }
-
+                                      
                                       for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) {
                                           // Only do work if the row group is not empty.
                                           if (*rowGroupingIt != *(rowGroupingIt + 1)) {
@@ -633,7 +634,7 @@ namespace storm {
                                               }
                                               
                                               for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) {
-                                                  if (filter(*sourceIt, *targetIt)) {
+                                                  if (f(*sourceIt, *targetIt)) {
                                                       *targetIt = *sourceIt;
                                                       if (choices != nullptr) {
                                                           *choiceIt = localChoice;
@@ -647,7 +648,7 @@ namespace storm {
                                           } else {
                                               // Compensate for the 'wrong' move forward in the loop header.
                                               --targetIt;
-
+                                              
                                               // Record dummy choice.
                                               if (choices != nullptr) {
                                                   *choiceIt = 0;
@@ -678,7 +679,7 @@ namespace storm {
                             *choiceIt = 0;
                         }
                         for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) {
-                            if (filter(*sourceIt, *targetIt)) {
+                            if (f(*sourceIt, *targetIt)) {
                                 *targetIt = *sourceIt;
                                 if (choices != nullptr) {
                                     *choiceIt = localChoice;
@@ -702,6 +703,8 @@ namespace storm {
 #endif
             }
             
+            
+            
             /*!
              * Reduces the given source vector by selecting the smallest element out of each row group.
              *
@@ -712,7 +715,7 @@ namespace storm {
              */
             template<class T>
             void reduceVectorMin(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
-                reduceVector<T>(source, target, rowGrouping, std::less<T>(), choices);
+                reduceVector<T, std::less<T>>(source, target, rowGrouping, choices);
             }
             
             /*!
@@ -725,7 +728,7 @@ namespace storm {
              */
             template<class T>
             void reduceVectorMax(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
-                reduceVector<T>(source, target, rowGrouping, std::greater<T>(), choices);
+                reduceVector<T, std::greater<T>>(source, target, rowGrouping, choices);
             }
             
             /*!
@@ -746,7 +749,6 @@ namespace storm {
                 }
             }
             
-            
             /*!
              * Compares the given elements and determines whether they are equal modulo the given precision. The provided flag
              * additionaly specifies whether the error is computed in relative or absolute terms.