From 3c844a487f6fe71761969b5ce186938288096c92 Mon Sep 17 00:00:00 2001 From: dehnert Date: Thu, 7 Sep 2017 22:15:10 +0200 Subject: [PATCH] some more optimizations --- .../prctl/helper/SparseMdpPrctlHelper.cpp | 20 +++++-- .../models/sparse/StandardRewardModel.cpp | 8 +-- src/storm/utility/vector.h | 52 ++++++++++--------- 3 files changed, 47 insertions(+), 33 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index af195ac4a..5e70c39e9 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -960,6 +960,8 @@ namespace storm { template std::unique_ptr SparseMdpPrctlHelper::computeConditionalProbabilities(OptimizationDirection dir, storm::storage::sparse::state_type initialState, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, storm::solver::MinMaxLinearEquationSolverFactory 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 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(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 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(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()); ++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(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 newTransitionMatrix = builder.build(); STORM_LOG_DEBUG("Transformed model has " << newTransitionMatrix.getRowGroupCount() << " states and " << newTransitionMatrix.getNonzeroEntryCount() << " transitions."); storm::storage::SparseMatrix newBackwardTransitions = newTransitionMatrix.transpose(true); + std::chrono::high_resolution_clock::time_point conditionalStart = std::chrono::high_resolution_clock::now(); std::vector 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(conditionalEnd - conditionalStart).count() << "ms."); return std::unique_ptr(new ExplicitQuantitativeCheckResult(initialState, dir == OptimizationDirection::Maximize ? goalProbabilities[numberOfStatesBeforeRelevantStates[initialState]] : storm::utility::one() - goalProbabilities[numberOfStatesBeforeRelevantStates[initialState]])); } diff --git a/src/storm/models/sparse/StandardRewardModel.cpp b/src/storm/models/sparse/StandardRewardModel.cpp index 4385c46e1..ccf9b5d71 100644 --- a/src/storm/models/sparse/StandardRewardModel.cpp +++ b/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(this->getStateActionRewardVector(), *weights, this->getStateRewardVector(), + storm::utility::vector::applyPointwiseTernary(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(this->optionalStateRewardVector.get(), *weights, this->optionalStateRewardVector.get(), [] (ValueType const& r, MatrixValueType const& w) { return w * r; } ); + storm::utility::vector::applyPointwise>(this->optionalStateRewardVector.get(), *weights, this->optionalStateRewardVector.get()); } } else { if (this->hasStateRewards()) { @@ -216,11 +216,11 @@ namespace storm { std::vector result; if (this->hasTransitionRewards()) { result = transitionMatrix.getPointwiseProductRowSumVector(this->getTransitionRewardMatrix()); - storm::utility::vector::applyPointwise(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * (resultElement + rewardElement); } ); + storm::utility::vector::applyPointwiseTernary(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * (resultElement + rewardElement); } ); } else { result = std::vector(transitionMatrix.getRowCount()); if (this->hasStateActionRewards()) { - storm::utility::vector::applyPointwise(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement, ValueType const& resultElement) { return weight * rewardElement; } ); + storm::utility::vector::applyPointwise(weights, this->getStateActionRewardVector(), result, [] (MatrixValueType const& weight, ValueType const& rewardElement) { return weight * rewardElement; } ); } } if (this->hasStateRewards()) { diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index a9ef8d0b0..89fbdf1bf 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -307,8 +307,8 @@ namespace storm { * @param secondOperand The second operand. * @param target The target vector. */ - template - void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, std::function const& function) { + template + void applyPointwiseTernary(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { #ifdef STORM_HAVE_INTELTBB tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range 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 - void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, std::function const& function) { + template + void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { #ifdef STORM_HAVE_INTELTBB tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range 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 - void applyPointwise(std::vector const& operand, std::vector& target, std::function const& function) { + template + void applyPointwise(std::vector const& operand, std::vector& target, Operation f = Operation()) { #ifdef STORM_HAVE_INTELTBB tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range 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 void addVectors(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target) { - applyPointwise(firstOperand, secondOperand, target, std::plus<>()); + applyPointwise>(firstOperand, secondOperand, target); } /*! @@ -397,7 +397,7 @@ namespace storm { */ template void subtractVectors(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target) { - applyPointwise(firstOperand, secondOperand, target, std::minus<>()); + applyPointwise>(firstOperand, secondOperand, target); } /*! @@ -409,7 +409,7 @@ namespace storm { */ template void multiplyVectorsPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target) { - applyPointwise(firstOperand, secondOperand, target, std::multiplies<>()); + applyPointwise>(firstOperand, secondOperand, target); } /*! @@ -421,7 +421,7 @@ namespace storm { */ template void divideVectorsPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target) { - applyPointwise(firstOperand, secondOperand, target, std::divides<>()); + applyPointwise>(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 - void reduceVector(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::function filter, std::vector* choices) { + template + void reduceVector(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices) { + Filter f; #ifdef STORM_HAVE_INTELTBB tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range 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 void reduceVectorMin(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) { - reduceVector(source, target, rowGrouping, std::less(), choices); + reduceVector>(source, target, rowGrouping, choices); } /*! @@ -725,7 +728,7 @@ namespace storm { */ template void reduceVectorMax(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) { - reduceVector(source, target, rowGrouping, std::greater(), choices); + reduceVector>(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.