diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ae07a267..978ec8b33 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -721,7 +721,12 @@ if (MSVC) target_link_libraries(storm "Dbghelp.lib") endif(MSVC) -# Print Cotire Usage Status + +############################################################# +## +## Cotire +## +############################################################# message (STATUS "StoRM - Using Cotire: ${STORM_USE_COTIRE}") if (STORM_USE_COTIRE) @@ -735,12 +740,23 @@ if (STORM_USE_COTIRE) #cotire(storm-performance-tests) endif() +############################################################# +## +## libc++abi +## +############################################################# # Link against libc++abi if requested. May be needed to build on Linux systems using clang. if (LINK_LIBCXXABI) message (STATUS "StoRM - Linking against libc++abi.") target_link_libraries(storm "c++abi") endif(LINK_LIBCXXABI) + +############################################################# +## +## Doxygen +## +############################################################# # Add a target to generate API documentation with Doxygen if(DOXYGEN_FOUND) set(CMAKE_DOXYGEN_OUTPUT_DIR "${CMAKE_CURRENT_BINARY_DIR}/doc") @@ -751,6 +767,11 @@ if(DOXYGEN_FOUND) add_custom_target(doc ${DOXYGEN_EXECUTABLE} "${CMAKE_CURRENT_BINARY_DIR}/Doxyfile" DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/Doxyfile" COMMENT "Generating API documentation with Doxygen" VERBATIM) endif(DOXYGEN_FOUND) +############################################################# +## +## memcheck targets +## +############################################################# add_custom_target(memcheck valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm -v --fix-deadlocks ${PROJECT_SOURCE_DIR}/examples/dtmc/crowds/crowds5_5.tra examples/dtmc/crowds/crowds5_5.lab DEPENDS storm) add_custom_target(memcheck-functional-tests valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm-functional-tests -v --fix-deadlocks DEPENDS storm-functional-tests) add_custom_target(memcheck-performance-tests valgrind --leak-check=full --show-reachable=yes ${PROJECT_BINARY_DIR}/storm-performance-tests -v --fix-deadlocks DEPENDS storm-performance-tests) diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp index e371ad8ce..498a2fc02 100644 --- a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp @@ -141,7 +141,7 @@ namespace storm { STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); // Compute the uniformized matrix. - storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates, statesWithProbabilityGreater0NonPsi, uniformizationRate); + storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates, statesWithProbabilityGreater0NonPsi, uniformizationRate); // Compute the vector that is to be added as a compensation for removing the absorbing states. storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * this->getModel().getTransitionMatrix() * psiStates.swapVariables(this->getModel().getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(this->getModel().getColumnVariables()) / this->getModel().getManager().getConstant(uniformizationRate); @@ -354,7 +354,24 @@ namespace storm { // Finally, compute the transient probabilities. std::vector<ValueType> result = SparseCtmcCslModelChecker<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, *this->linearEquationSolverFactory); - return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), odd, result)); + return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result))); + } + + template<storm::dd::DdType DdType, class ValueType> + std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula); + SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + + storm::dd::Add<DdType> probabilityMatrix = this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + + // Create ODD for the translation. + storm::dd::Odd<DdType> odd(this->getModel().getReachableStates()); + + storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd); + std::vector<ValueType> explicitExitRateVector = this->getModel().getExitRateVector().template toVector<ValueType>(odd); + + std::vector<ValueType> result = SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), &explicitExitRateVector, qualitative, *this->linearEquationSolverFactory); + return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result))); } // Explicitly instantiate the model checker. diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.h b/src/modelchecker/csl/HybridCtmcCslModelChecker.h index f0e619533..b2a94d51e 100644 --- a/src/modelchecker/csl/HybridCtmcCslModelChecker.h +++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.h @@ -23,7 +23,8 @@ namespace storm { virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; - + virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + protected: storm::models::symbolic::Ctmc<DdType> const& getModel() const override; diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index ecd550f4c..9b6818884 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -12,6 +12,8 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/storage/StronglyConnectedComponentDecomposition.h" + #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/NotImplementedException.h" @@ -189,7 +191,7 @@ namespace storm { std::vector<ValueType> newSubresult = std::vector<ValueType>(relevantStates.getNumberOfSetBits()); storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one<ValueType>()); - + // Then compute the transient probabilities of being in such a state after t time units. For this, // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. uniformizationRate = storm::utility::zero<ValueType>(); @@ -344,7 +346,7 @@ namespace storm { weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; storm::utility::vector::applyPointwise(result, values, result, addAndScale); } - + return result; } @@ -406,7 +408,7 @@ namespace storm { } uniformizationRate *= 1.02; STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - + storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector()); // Compute the total state reward vector. @@ -436,6 +438,22 @@ namespace storm { return result; } + template<class ValueType> + storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<ValueType>::computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) { + storm::storage::SparseMatrix<ValueType> generatorMatrix(rateMatrix, true); + + // Place the negative exit rate on the diagonal. + for (uint_fast64_t row = 0; row < generatorMatrix.getRowCount(); ++row) { + for (auto& entry : generatorMatrix.getRow(row)) { + if (entry.getColumn() == row) { + entry.setValue(-exitRates[row]); + } + } + } + + return generatorMatrix; + } + template<class ValueType> std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula()); @@ -455,6 +473,240 @@ namespace storm { return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative))); } + template<class ValueType> + std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + + storm::storage::SparseMatrix<ValueType> probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory))); + } + + template<typename ValueType> + std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) { + // If there are no goal states, we avoid the computation and directly return zero. + uint_fast64_t numOfStates = transitionMatrix.getRowCount(); + if (psiStates.empty()) { + return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>()); + } + + // Likewise, if all bits are set, we can avoid the computation. + if (psiStates.full()) { + return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>()); + } + + // Start by decomposing the DTMC into its BSCCs. + storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowCount(), true), false, true); + + STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs."); + + // Get some data members for convenience. + ValueType one = storm::utility::one<ValueType>(); + ValueType zero = storm::utility::zero<ValueType>(); + + // Prepare the vector holding the LRA values for each of the BSCCs. + std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero); + + // First we check which states are in BSCCs. + storm::storage::BitVector statesInBsccs(numOfStates); + storm::storage::BitVector firstStatesInBsccs(numOfStates); + + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + + // Gather information for later use. + bool first = true; + for (auto const& state : bscc) { + statesInBsccs.set(state); + if (first) { + firstStatesInBsccs.set(state); + } + first = false; + } + } + storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; + + STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs."); + + // Prepare a vector holding the index within all states that are in BSCCs for every state. + std::vector<uint_fast64_t> indexInStatesInBsccs; + + // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC. + std::vector<uint_fast64_t> stateToBsccIndexMap; + + if (!statesInBsccs.empty()) { + firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; + + // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. + storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); + + // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a + // multiplication from the right by transposing the system. + bsccEquationSystem = bsccEquationSystem.transpose(false, true); + + // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states. + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + indexInStatesInBsccs.reserve(transitionMatrix.getRowCount()); + for (auto index : statesInBsccs) { + while (lastIndex <= index) { + indexInStatesInBsccs.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + + stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits()); + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + for (auto const& state : bscc) { + stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex; + } + } + + // Now build the final equation system matrix, the initial guess and the right-hand side in one go. + std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); + storm::storage::SparseMatrixBuilder<ValueType> builder; + for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { + + // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the + // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the + // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal. + if (firstStatesInBsccs.get(row)) { + uint_fast64_t requiredBscc = stateToBsccIndexMap[row]; + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc]; + + for (auto const& state : bscc) { + builder.addNextValue(row, indexInStatesInBsccs[state], one); + } + + bsccEquationSystemRightSide[row] = one; + + } else { + // Otherwise, we copy the row, and subtract 1 from the diagonal. + for (auto& entry : bsccEquationSystem.getRow(row)) { + if (entry.getColumn() == row) { + builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); + } else { + builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } + } + } + + } + + // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC. + std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); + } + } + + bsccEquationSystem = builder.build(); + + { + std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem); + solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); + } + + // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. + if (exitRateVector != nullptr) { + std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero); + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { + bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]); + } + + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { + bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]]; + } + } + // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + if (psiStates.get(state)) { + bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += bsccEquationSystemSolution[indexInStatesInBsccs[state]]; + } + } + } + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); + } + } else { + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values + // directly (based on whether or not the contained state is a psi state). + if (psiStates.get(*bscc.begin())) { + bsccLra[bsccIndex] = 1; + } + } + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); + } + } + + std::vector<ValueType> rewardSolution; + if (!statesNotInBsccs.empty()) { + // Calculate LRA for states not in bsccs as expected reachability rewards. + // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a + // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability + // of the BSCC. + + std::vector<ValueType> rewardRightSide; + rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); + + for (auto state : statesNotInBsccs) { + ValueType reward = zero; + for (auto entry : transitionMatrix.getRow(state)) { + if (statesInBsccs.get(entry.getColumn())) { + reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]]; + } + } + rewardRightSide.push_back(reward); + } + + storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); + rewardEquationSystemMatrix.convertToEquationSystem(); + + rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one); + + { + std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix); + solver->solveEquationSystem(rewardSolution, rewardRightSide); + } + } + + // Fill the result vector. + std::vector<ValueType> result(numOfStates); + auto rewardSolutionIter = rewardSolution.begin(); + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + result[state] = bsccLra[bsccIndex]; + } + } + for (auto state : statesNotInBsccs) { + STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); + // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th + // entry in rewardSolution we can just take the next value from the iterator. + result[state] = *rewardSolutionIter; + ++rewardSolutionIter; + } + + return result; + } + + // Explicitly instantiate the model checker. template class SparseCtmcCslModelChecker<double>; diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h index 5108a183f..af1c9394a 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -13,10 +13,14 @@ namespace storm { template<storm::dd::DdType DdType, typename ValueType> class HybridCtmcCslModelChecker; + template<typename ValueType> + class SparseDtmcPrctlModelChecker; + template<class ValueType> class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> { public: friend class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, ValueType>; + friend class SparseDtmcPrctlModelChecker<ValueType>; explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model); explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory); @@ -29,6 +33,7 @@ namespace storm { virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; protected: storm::models::sparse::Ctmc<ValueType> const& getModel() const override; @@ -39,7 +44,7 @@ namespace storm { static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory); std::vector<ValueType> computeInstantaneousRewardsHelper(double timeBound) const; std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const; - + static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory); /*! * Computes the matrix representing the transitions of the uniformized CTMC. * @@ -73,8 +78,18 @@ namespace storm { * * @param rateMatrix The rate matrix. * @param exitRates The exit rate vector. + * @return The †ransition matrix of the embedded DTMC. */ static storm::storage::SparseMatrix<ValueType> computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates); + + /*! + * Converts the given rate-matrix into the generator matrix. + * + * @param rateMatrix The rate matrix. + * @param exitRates The exit rate vector. + * @return The generator matrix. + */ + static storm::storage::SparseMatrix<ValueType> computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates); // An object that is used for solving linear equations and performing matrix-vector multiplication. std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory; diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp index 25a02adff..3b61ded53 100644 --- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp @@ -1,4 +1,5 @@ #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h" +#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "src/storage/dd/CuddOdd.h" @@ -301,6 +302,20 @@ namespace storm { return this->template getModelAs<storm::models::symbolic::Dtmc<DdType>>(); } + template<storm::dd::DdType DdType, class ValueType> + std::unique_ptr<CheckResult> HybridDtmcPrctlModelChecker<DdType, ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula); + SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + + // Create ODD for the translation. + storm::dd::Odd<DdType> odd(this->getModel().getReachableStates()); + + storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = this->getModel().getTransitionMatrix().toMatrix(odd, odd); + + std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), qualitative, *this->linearEquationSolverFactory); + return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result))); + } + template class HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double>; } } \ No newline at end of file diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h index b3ab51279..8964d29d7 100644 --- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h @@ -26,6 +26,7 @@ namespace storm { virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; protected: storm::models::symbolic::Dtmc<DdType> const& getModel() const override; diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index f16951d12..d49075faa 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -8,11 +8,11 @@ #include "src/utility/graph.h" #include "src/utility/solver.h" +#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "src/exceptions/InvalidStateException.h" -#include "src/storage/StronglyConnectedComponentDecomposition.h" #include "src/exceptions/InvalidPropertyException.h" @@ -302,232 +302,18 @@ namespace storm { ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative))); } - - template<typename ValueType> - std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const { - // If there are no goal states, we avoid the computation and directly return zero. - auto numOfStates = this->getModel().getNumberOfStates(); - if (psiStates.empty()) { - return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>()); - } - - // Likewise, if all bits are set, we can avoid the computation and set. - if ((~psiStates).empty()) { - return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>()); - } - - // Start by decomposing the DTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(this->getModel(), false, true); - - // Get some data members for convenience. - typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix(); - ValueType one = storm::utility::one<ValueType>(); - ValueType zero = storm::utility::zero<ValueType>(); - - // First we check which states are in BSCCs. We use this later to speed up reachability analysis - storm::storage::BitVector statesInBsccs(numOfStates); - - std::vector<uint_fast64_t> stateToBsccIndexMap(transitionMatrix.getColumnCount()); - - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // Gather information for later use. - for (auto const& state : bscc) { - statesInBsccs.set(state); - stateToBsccIndexMap[state] = currentBsccIndex; - } - } - - storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; - - // calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of the transposed transition matrix for the bsccs - storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); - - //subtract identity matrix - for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { - for (auto& entry : bsccEquationSystem.getRow(row)) { - if (entry.getColumn() == row) { - entry.setValue(entry.getValue() - one); - } - } - } - //now transpose, this internally removes all explicit zeros from the matrix that where introduced when subtracting the identity matrix - bsccEquationSystem = bsccEquationSystem.transpose(); - - std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); - std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); - { - auto solver = this->linearEquationSolverFactory->create(bsccEquationSystem); - solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); - } - - //calculate LRA Value for each BSCC from steady state distribution in BSCCs - // we have to do some scaling, as the probabilities for each BSCC have to sum up to one, which they don't necessarily do in the solution of the equation system - std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero); - std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero); - size_t i = 0; - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { - if (psiStates.get(*stateIter)) { - bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; - } - bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; - } - for (i = 0; i < bsccLra.size(); ++i) { - bsccLra[i] /= bsccTotalValue[i]; - } - - //calculate LRA for states not in bsccs as expected reachability rewards - //target states are states in bsccs, transition reward is the lra of the bscc for each transition into a bscc and 0 otherwise - //this corresponds to sum of LRAs in BSCC weighted by the reachability probability of the BSCC - - std::vector<ValueType> rewardRightSide; - rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); - - for (auto state : statesNotInBsccs) { - ValueType reward = zero; - for (auto entry : transitionMatrix.getRow(state)) { - if (statesInBsccs.get(entry.getColumn())) { - reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; - } - } - rewardRightSide.push_back(reward); - } - - storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); - rewardEquationSystemMatrix.convertToEquationSystem(); - - std::vector<ValueType> rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one); - - { - auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix); - solver->solveEquationSystem(rewardSolution, rewardRightSide); - } - - // now fill the result vector - std::vector<ValueType> result(numOfStates); - - auto rewardSolutionIter = rewardSolution.begin(); - for (size_t state = 0; state < numOfStates; ++state) { - if (statesInBsccs.get(state)) { - //assign the value of the bscc the state is in - result[state] = bsccLra[stateToBsccIndexMap[state]]; - } else { - assert(rewardSolutionIter != rewardSolution.end()); - //take the value from the reward computation - //since the n-th state not in any bscc is the n-th entry in rewardSolution we can just take the next value from the iterator - result[state] = *rewardSolutionIter; - ++rewardSolutionIter; - } - } - - return result; - - //old implementeation - - //now we calculate the long run average for each BSCC in isolation and compute the weighted contribution of the BSCC to the LRA value of all states - //for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - // storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // storm::storage::BitVector statesInThisBscc(numOfStates); - // for (auto const& state : bscc) { - // statesInThisBscc.set(state); - // } - - // //ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); - // ValueType lraForThisBscc = bsccLra[currentBsccIndex]; - - // //the LRA value of a BSCC contributes to the LRA value of a state with the probability of reaching that BSCC from that state - // //thus we add Prob(Eventually statesInThisBscc) * lraForThisBscc to the result vector - - // //the reachability probabilities will be zero in other BSCCs, thus we can set the left operand of the until formula to statesNotInBsccs as an optimization - // std::vector<ValueType> reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false); - // - // storm::utility::vector::scaleVectorInPlace<ValueType>(reachProbThisBscc, lraForThisBscc); - // storm::utility::vector::addVectorsInPlace<ValueType>(result, reachProbThisBscc); - //} - - //return result; - } - - template<typename ValueType> - std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { - std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula); - ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - - return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative))); - } - + template<typename ValueType> + std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + + return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory))); + } + template<typename ValueType> - ValueType SparseDtmcPrctlModelChecker<ValueType>::computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) { - //if size is one we can immediately derive the result - if (bscc.size() == 1){ - if (psiStates.get(*(bscc.begin()))) { - return storm::utility::one<ValueType>(); - } else{ - return storm::utility::zero<ValueType>(); - } - } - - storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount()); - subsystem.set(bscc.begin(), bscc.end()); - - //we now have to solve ((P')^t - I) * x = 0, where P' is the submatrix of transitionMatrix, - // ^t means transose, and I is the identity matrix. - - storm::storage::SparseMatrix<ValueType> subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem, true); - subsystemMatrix = subsystemMatrix.transpose(); - - // subtract 1 on the diagonal and at the same time add a row with all ones to enforce that the result of the equation system is unique - storm::storage::SparseMatrixBuilder<ValueType> equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount()); - ValueType one = storm::utility::one<ValueType>(); - ValueType zero = storm::utility::zero<ValueType>(); - bool foundDiagonalElement = false; - for (uint_fast64_t row = 0; row < subsystemMatrix.getRowCount(); ++row) { - for (auto& entry : subsystemMatrix.getRowGroup(row)) { - if (entry.getColumn() == row) { - equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue() - one); - foundDiagonalElement = true; - } else { - equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue()); - } - } - - // Throw an exception if a row did not have an element on the diagonal. - STORM_LOG_THROW(foundDiagonalElement, storm::exceptions::InvalidOperationException, "Internal Error, no diagonal entry found."); - } - for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) { - equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one); - } - equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero); - subsystemMatrix = equationSystemBuilder.build(); - - // create x and b for the equation system. setting the last entry of b to one enforces the sum over the unique solution vector is one - std::vector<ValueType> steadyStateDist(subsystemMatrix.getRowCount(), zero); - std::vector<ValueType> b(subsystemMatrix.getRowCount(), zero); - b[subsystemMatrix.getRowCount() - 1] = one; - - - { - auto solver = linearEquationSolverFactory.create(subsystemMatrix); - solver->solveEquationSystem(steadyStateDist, b); - } - - //remove the last entry of the vector as it was just there to enforce the unique solution - steadyStateDist.pop_back(); - - //calculate the fraction we spend in psi states on the long run - std::vector<ValueType> steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero); - storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist); - - ValueType result = zero; - - for (auto value : steadyStateForPsiStates) { - result += value; - } - - return result; + std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) { + return SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory); } template<typename ValueType> diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index b55f00faa..d37023948 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -9,12 +9,16 @@ namespace storm { namespace modelchecker { + template<storm::dd::DdType DdType, typename ValueType> + class HybridDtmcPrctlModelChecker; + // Forward-declare CTMC model checker so we can make it a friend. template<typename ValueType> class SparseCtmcCslModelChecker; template<class ValueType> class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker<ValueType> { public: + friend class HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, ValueType>; friend class SparseCtmcCslModelChecker<ValueType>; explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model); @@ -41,9 +45,7 @@ namespace storm { std::vector<ValueType> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const; std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const; static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative); - std::vector<ValueType> computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const; - - static ValueType computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory); + static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory); // An object that is used for retrieving linear equation solvers. std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory; diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp new file mode 100644 index 000000000..022eca0af --- /dev/null +++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp @@ -0,0 +1,263 @@ +#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h" + +#include "src/storage/dd/CuddOdd.h" + +#include "src/utility/macros.h" +#include "src/utility/graph.h" + +#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h" +#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h" +#include "src/modelchecker/results/HybridQuantitativeCheckResult.h" + +#include "src/exceptions/InvalidStateException.h" +#include "src/exceptions/InvalidPropertyException.h" + +namespace storm { + namespace modelchecker { + template<storm::dd::DdType DdType, typename ValueType> + SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) { + // Intentionally left empty. + } + + template<storm::dd::DdType DdType, typename ValueType> + SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>()) { + // Intentionally left empty. + } + + template<storm::dd::DdType DdType, typename ValueType> + bool SymbolicDtmcPrctlModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const { + return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula(); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) { + // We need to identify the states which have to be taken out of the matrix, i.e. all states that have + // probability 0 and 1 of satisfying the until-formula. + std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates); + storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates(); + + // Perform some logging. + STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states."); + STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states."); + STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); + + // Check whether we need to compute exact probabilities for some states. + if (qualitative) { + // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1. + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5))); + } else { + // If there are maybe states, we need to solve an equation system. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd<DdType> odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd(); + prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn; + subvector = subvector.sumAbstract(model.getColumnVariables()); + + // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed + // for solving the equation system (i.e. compute (I-A)). + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; + + // Solve the equation system. + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector); + + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + result)); + } else { + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd())); + } + } + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula()); + std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula()); + SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + return this->computeUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula()); + SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector()))); + } + + template<storm::dd::DdType DdType, typename ValueType> + storm::dd::Add<DdType> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) { + storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd(); + return result.sumAbstract(model.getColumnVariables()); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); + std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula()); + std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula()); + SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + return this->computeBoundedUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) { + // We need to identify the states which have to be taken out of the matrix, i.e. all states that have + // probability 0 or 1 of satisfying the until-formula. + storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound); + storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates(); + + // If there are maybe states, we need to perform matrix-vector multiplications. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd<DdType> odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables()); + + // Finally cut away all columns targeting non-maybe states. + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + + // Perform the matrix-vector multiplication. + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &subvector, stepBound); + + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd() + result)); + } else { + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd())); + } + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); + return this->computeCumulativeRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) { + // Only compute the result if the model has at least one reward this->getModel(). + STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Compute the reward vector to add in each step based on the available reward models. + storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero(); + if (model.hasTransitionRewards()) { + totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables()); + } + + // Perform the matrix-vector multiplication. + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &totalRewardVector, stepBound); + + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result)); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); + return this->computeInstantaneousRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) { + // Only compute the result if the model has at least one reward this->getModel(). + STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Perform the matrix-vector multiplication. + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getStateRewardVector(), nullptr, stepBound); + + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result)); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) { + std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula()); + SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>(); + return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative); + } + + template<storm::dd::DdType DdType, typename ValueType> + std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative) { + + // Only compute the result if there is at least one reward model. + STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Determine which states have a reward of infinity by definition. + storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates); + infinityStates = !infinityStates && model.getReachableStates(); + storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates(); + STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states."); + STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states."); + STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); + + // Check whether we need to compute exact rewards for some states. + if (qualitative) { + // Set the values for all maybe-states to 1 to indicate that their reward values + // are neither 0 nor infinity. + return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>()))); + } else { + // If there are maybe states, we need to solve an equation system. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd<DdType> odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the state reward vector to use in the computation. + storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero(); + if (transitionRewardMatrix) { + subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables()); + } + + // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed + // for solving the equation system (i.e. compute (I-A)). + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; + + // Solve the equation system. + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector); + + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + result)); + } else { + return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()))); + } + } + } + + template<storm::dd::DdType DdType, typename ValueType> + storm::models::symbolic::Dtmc<DdType> const& SymbolicDtmcPrctlModelChecker<DdType, ValueType>::getModel() const { + return this->template getModelAs<storm::models::symbolic::Dtmc<DdType>>(); + } + + template class SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double>; + } +} \ No newline at end of file diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h new file mode 100644 index 000000000..9dca0335b --- /dev/null +++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h @@ -0,0 +1,49 @@ +#ifndef STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_ +#define STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_ + +#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h" +#include "src/models/symbolic/Dtmc.h" +#include "src/utility/solver.h" + +namespace storm { + namespace modelchecker { + template<storm::dd::DdType DdType, typename ValueType> + class SymbolicCtmcCslModelChecker; + + template<storm::dd::DdType DdType, typename ValueType> + class SymbolicDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> { + public: + friend class SymbolicCtmcCslModelChecker<DdType, ValueType>; + + explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model); + explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory); + + // The implemented methods of the AbstractModelChecker interface. + virtual bool canHandle(storm::logic::Formula const& formula) const override; + virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override; + + protected: + storm::models::symbolic::Dtmc<DdType> const& getModel() const override; + + private: + // The methods that perform the actual checking. + static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory); + static storm::dd::Add<DdType> computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates); + static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory); + static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory); + static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory); + static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative); + + // An object that is used for retrieving linear equation solvers. + std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>> linearEquationSolverFactory; + }; + + } // namespace modelchecker +} // namespace storm + +#endif /* STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_ */ diff --git a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp index 4fa0dbafe..1ddd761d6 100644 --- a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp +++ b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp @@ -36,6 +36,11 @@ namespace storm { return std::unique_ptr<CheckResult>(new SymbolicQualitativeCheckResult<Type>(model.getReachableStates(), model.getStates(stateFormula.getLabel()))); } + template<storm::dd::DdType Type> + std::unique_ptr<CheckResult> SymbolicPropositionalModelChecker<Type>::checkAtomicExpressionFormula(storm::logic::AtomicExpressionFormula const& stateFormula) { + return std::unique_ptr<CheckResult>(new SymbolicQualitativeCheckResult<Type>(model.getReachableStates(), model.getStates(stateFormula.getExpression()))); + } + template<storm::dd::DdType Type> storm::models::symbolic::Model<Type> const& SymbolicPropositionalModelChecker<Type>::getModel() const { return model; diff --git a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.h b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.h index 96dcea601..ec3433dec 100644 --- a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.h +++ b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.h @@ -16,7 +16,8 @@ namespace storm { virtual bool canHandle(storm::logic::Formula const& formula) const override; virtual std::unique_ptr<CheckResult> checkBooleanLiteralFormula(storm::logic::BooleanLiteralFormula const& stateFormula) override; virtual std::unique_ptr<CheckResult> checkAtomicLabelFormula(storm::logic::AtomicLabelFormula const& stateFormula) override; - + virtual std::unique_ptr<CheckResult> checkAtomicExpressionFormula(storm::logic::AtomicExpressionFormula const& stateFormula) override; + protected: /*! * Retrieves the model associated with this model checker instance. diff --git a/src/modelchecker/reachability/CollectConstraints.h b/src/modelchecker/reachability/CollectConstraints.h deleted file mode 100644 index d77273770..000000000 --- a/src/modelchecker/reachability/CollectConstraints.h +++ /dev/null @@ -1,63 +0,0 @@ -/** - * @file: CollectConstraints.h - * @author: Sebastian Junges - * - * @since October 8, 2014 - */ - -#pragma once -#include "src/models/Dtmc.h" - -namespace storm { -namespace modelchecker { - namespace reachability { - template<typename ValueType> - class CollectConstraints - { - private: - std::unordered_set<carl::Constraint<ValueType>> wellformedConstraintSet; - std::unordered_set<carl::Constraint<ValueType>> graphPreservingConstraintSet; - storm::utility::ConstantsComparator<ValueType> comparator; - - public: - std::unordered_set<carl::Constraint<ValueType>> const& wellformedConstraints() const { - return this->wellformedConstraintSet; - } - - std::unordered_set<carl::Constraint<ValueType>> const& graphPreservingConstraints() const { - return this->graphPreservingConstraintSet; - } - - void process(storm::models::Dtmc<ValueType> const& dtmc) - { - for(uint_fast64_t state = 0; state < dtmc.getNumberOfStates(); ++state) - { - ValueType sum; - assert(comparator.isZero(sum)); - for(auto const& transition : dtmc.getRows(state)) - { - sum += transition.getValue(); - if(!transition.getValue().isConstant()) - { - wellformedConstraintSet.emplace(transition.getValue() - 1, storm::CompareRelation::LEQ); - wellformedConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GEQ); - graphPreservingConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GT); - } - } - assert(!comparator.isConstant(sum) || comparator.isOne(sum)); - if(!sum.isConstant()) { - wellformedConstraintSet.emplace(sum - 1, storm::CompareRelation::EQ); - } - - } - } - - void operator()(storm::models::Dtmc<ValueType> const& dtmc) - { - process(dtmc); - } - - }; - } -} -} diff --git a/src/modelchecker/reachability/DirectEncoding.h b/src/modelchecker/reachability/DirectEncoding.h deleted file mode 100644 index 69ddeeb28..000000000 --- a/src/modelchecker/reachability/DirectEncoding.h +++ /dev/null @@ -1,80 +0,0 @@ -/** - * @file: DirectEncoding.h - * @author: Sebastian Junges - * - * @since April 8, 2014 - */ - -#pragma once - -#ifdef STORM_HAVE_CARL -#include <carl/io/WriteTosmt2Stream.h> - -namespace storm -{ - namespace modelchecker - { - namespace reachability - { - class DirectEncoding - { - public: - template<typename T> - std::string encodeAsSmt2(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& oneStepProbabilities, std::set<carl::Variable> const& parameters, storm::storage::BitVector const& initialStates, typename T::CoeffType const& threshold, bool lessequal = false) { - - carl::io::WriteTosmt2Stream smt2; - uint_fast64_t numberOfStates = transitionMatrix.getRowCount(); - carl::VariablePool& vpool = carl::VariablePool::getInstance(); - std::vector<carl::Variable> stateVars; - for (carl::Variable const& p : parameters) { - smt2 << ("parameter_bound_" + vpool.getName(p)); - smt2 << carl::io::smt2node::AND; - smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(p), carl::CompareRelation::GT); - smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(p) - Polynomial::PolyType(1), carl::CompareRelation::LT); - smt2 << carl::io::smt2node::CLOSENODE; - } - - for (uint_fast64_t state = 0; state < numberOfStates; ++state) { - carl::Variable stateVar = vpool.getFreshVariable("s_" + std::to_string(state)); - stateVars.push_back(stateVar); - smt2 << ("state_bound_" + std::to_string(state)); - smt2 << carl::io::smt2node::AND; - smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar), carl::CompareRelation::GT); - smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar) - Polynomial::PolyType(1), carl::CompareRelation::LT); - smt2 << carl::io::smt2node::CLOSENODE; - } - - smt2.setAutomaticLineBreaks(true); - Polynomial::PolyType initStateReachSum; - for (uint_fast64_t state = 0; state < numberOfStates; ++state) { - T reachpropPol; - for (auto const& transition : transitionMatrix.getRow(state)) { -// reachpropPol += transition.getValue() * stateVars[transition.getColumn()]; - } - reachpropPol += oneStepProbabilities[state]; - smt2 << ("transition_" + std::to_string(state)); -// smt2 << carl::Constraint<Polynomial::PolyType>(reachpropPol - stateVars[state], carl::CompareRelation::EQ); - } - - smt2 << ("reachability"); - - carl::CompareRelation thresholdRelation = lessequal ? carl::CompareRelation::LEQ : carl::CompareRelation::GEQ; - smt2 << carl::io::smt2node::OR; - for (uint_fast64_t state : initialStates) { - smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVars[state]) - threshold, thresholdRelation); - } - smt2 << carl::io::smt2node::CLOSENODE; - - smt2 << carl::io::smt2flag::CHECKSAT; - smt2 << carl::io::smt2flag::MODEL; - smt2 << carl::io::smt2flag::UNSAT_CORE; - std::stringstream strm; - strm << smt2; - return strm.str(); - } - }; - } - } -} - -#endif \ No newline at end of file diff --git a/src/modelchecker/results/ExplicitQualitativeCheckResult.cpp b/src/modelchecker/results/ExplicitQualitativeCheckResult.cpp index 79bb34885..121721846 100644 --- a/src/modelchecker/results/ExplicitQualitativeCheckResult.cpp +++ b/src/modelchecker/results/ExplicitQualitativeCheckResult.cpp @@ -29,13 +29,19 @@ namespace storm { // Intentionally left empty. } - void ExplicitQualitativeCheckResult::performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, std::function<bool (bool, bool)> const& function) { + void ExplicitQualitativeCheckResult::performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, bool logicalAnd) { STORM_LOG_THROW(second.isExplicitQualitativeCheckResult(), storm::exceptions::InvalidOperationException, "Cannot perform logical 'and' on check results of incompatible type."); STORM_LOG_THROW(first.isResultForAllStates() == second.isResultForAllStates(), storm::exceptions::InvalidOperationException, "Cannot perform logical 'and' on check results of incompatible type."); ExplicitQualitativeCheckResult const& secondCheckResult = static_cast<ExplicitQualitativeCheckResult const&>(second); if (first.isResultForAllStates()) { - boost::get<vector_type>(first.truthValues) &= boost::get<vector_type>(secondCheckResult.truthValues); + if (logicalAnd) { + boost::get<vector_type>(first.truthValues) &= boost::get<vector_type>(secondCheckResult.truthValues); + } else { + boost::get<vector_type>(first.truthValues) |= boost::get<vector_type>(secondCheckResult.truthValues); + } } else { + std::function<bool (bool, bool)> function = logicalAnd ? [] (bool a, bool b) { return a && b; } : [] (bool a, bool b) { return a || b; }; + map_type& map1 = boost::get<map_type>(first.truthValues); map_type const& map2 = boost::get<map_type>(secondCheckResult.truthValues); for (auto& element1 : map1) { @@ -53,12 +59,12 @@ namespace storm { } QualitativeCheckResult& ExplicitQualitativeCheckResult::operator&=(QualitativeCheckResult const& other) { - performLogicalOperation(*this, other, [] (bool a, bool b) { return a && b; }); + performLogicalOperation(*this, other, true); return *this; } QualitativeCheckResult& ExplicitQualitativeCheckResult::operator|=(QualitativeCheckResult const& other) { - performLogicalOperation(*this, other, [] (bool a, bool b) { return a || b; }); + performLogicalOperation(*this, other, false); return *this; } diff --git a/src/modelchecker/results/ExplicitQualitativeCheckResult.h b/src/modelchecker/results/ExplicitQualitativeCheckResult.h index 8ab3e473d..abce327b4 100644 --- a/src/modelchecker/results/ExplicitQualitativeCheckResult.h +++ b/src/modelchecker/results/ExplicitQualitativeCheckResult.h @@ -51,7 +51,7 @@ namespace storm { virtual void filter(QualitativeCheckResult const& filter) override; private: - static void performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, std::function<bool (bool, bool)> const& function); + static void performLogicalOperation(ExplicitQualitativeCheckResult& first, QualitativeCheckResult const& second, bool logicalAnd); // The values of the quantitative check result. boost::variant<vector_type, map_type> truthValues; diff --git a/src/models/ModelType.cpp b/src/models/ModelType.cpp index 06edc9de2..5192a8aa7 100644 --- a/src/models/ModelType.cpp +++ b/src/models/ModelType.cpp @@ -8,6 +8,7 @@ namespace storm { case ModelType::Ctmc: os << "CTMC"; break; case ModelType::Mdp: os << "MDP"; break; case ModelType::MarkovAutomaton: os << "Markov Automaton"; break; + case ModelType::S2pg: os << "S2PG"; break; } return os; } diff --git a/src/models/ModelType.h b/src/models/ModelType.h index b67a03685..1c386030f 100644 --- a/src/models/ModelType.h +++ b/src/models/ModelType.h @@ -7,7 +7,7 @@ namespace storm { namespace models { // All supported model types. enum class ModelType { - Dtmc, Ctmc, Mdp, MarkovAutomaton + Dtmc, Ctmc, Mdp, MarkovAutomaton, S2pg }; std::ostream& operator<<(std::ostream& os, ModelType const& type); diff --git a/src/models/sparse/Ctmc.cpp b/src/models/sparse/Ctmc.cpp index e495221b1..5217396ce 100644 --- a/src/models/sparse/Ctmc.cpp +++ b/src/models/sparse/Ctmc.cpp @@ -10,7 +10,7 @@ namespace storm { Ctmc<ValueType>::Ctmc(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : DeterministicModel<ValueType>(storm::models::ModelType::Ctmc, rateMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { exitRates = createExitRateVector(this->getTransitionMatrix()); } @@ -19,7 +19,7 @@ namespace storm { Ctmc<ValueType>::Ctmc(storm::storage::SparseMatrix<ValueType>&& rateMatrix, storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : DeterministicModel<ValueType>(storm::models::ModelType::Ctmc, std::move(rateMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { // It is important to refer to the transition matrix here, because the given rate matrix has been move elsewhere. exitRates = createExitRateVector(this->getTransitionMatrix()); diff --git a/src/models/sparse/Ctmc.h b/src/models/sparse/Ctmc.h index 53be5048e..71e1eb074 100644 --- a/src/models/sparse/Ctmc.h +++ b/src/models/sparse/Ctmc.h @@ -29,7 +29,7 @@ namespace storm { Ctmc(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -43,7 +43,7 @@ namespace storm { Ctmc(storm::storage::SparseMatrix<ValueType>&& rateMatrix, storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); Ctmc(Ctmc<ValueType> const& ctmc) = default; Ctmc& operator=(Ctmc<ValueType> const& ctmc) = default; diff --git a/src/models/sparse/DeterministicModel.cpp b/src/models/sparse/DeterministicModel.cpp index d24353d02..68d18ba77 100644 --- a/src/models/sparse/DeterministicModel.cpp +++ b/src/models/sparse/DeterministicModel.cpp @@ -11,7 +11,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : Model<ValueType>(modelType, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { // Intentionally left empty. } @@ -22,7 +22,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : Model<ValueType>(modelType, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { // Intentionally left empty. } diff --git a/src/models/sparse/DeterministicModel.h b/src/models/sparse/DeterministicModel.h index 940048dc6..19ecb4a44 100644 --- a/src/models/sparse/DeterministicModel.h +++ b/src/models/sparse/DeterministicModel.h @@ -29,7 +29,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -46,7 +46,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); DeterministicModel(DeterministicModel const& other) = default; DeterministicModel& operator=(DeterministicModel const& other) = default; diff --git a/src/models/sparse/Dtmc.cpp b/src/models/sparse/Dtmc.cpp index d522743a3..f2f197428 100644 --- a/src/models/sparse/Dtmc.cpp +++ b/src/models/sparse/Dtmc.cpp @@ -12,7 +12,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : DeterministicModel<ValueType>(storm::models::ModelType::Dtmc, probabilityMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid."); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist."); @@ -22,7 +22,7 @@ namespace storm { Dtmc<ValueType>::Dtmc(storm::storage::SparseMatrix<ValueType>&& probabilityMatrix, storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : DeterministicModel<ValueType>(storm::models::ModelType::Dtmc, std::move(probabilityMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid."); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist."); @@ -41,7 +41,7 @@ namespace storm { // storm::models::sparse::StateLabeling(this->getStateLabeling(), subSysStates), // boost::optional<std::vector<ValueType>>(), // boost::optional<storm::storage::SparseMatrix<ValueType>>(), - // boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + // boost::optional<std::vector<LabelSet>>()); // } // // // Does the vector have the right size? @@ -164,16 +164,16 @@ namespace storm { // newTransitionRewards = newTransRewardsBuilder.build(); // } // - // boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> newChoiceLabels; + // boost::optional<std::vector<LabelSet>> newChoiceLabels; // if(this->hasChoiceLabeling()) { // // // Get the choice label sets and move the needed values to the front. - // std::vector<boost::container::flat_set<uint_fast64_t>> newChoice(this->getChoiceLabeling()); + // std::vector<LabelSet> newChoice(this->getChoiceLabeling()); // storm::utility::vector::selectVectorValues(newChoice, subSysStates, this->getChoiceLabeling()); // // // Throw away all values after the last state and set the choice label set for s_b as empty. // newChoice.resize(newStateCount); - // newChoice[newStateCount - 1] = boost::container::flat_set<uint_fast64_t>(); + // newChoice[newStateCount - 1] = LabelSet(); // // newChoiceLabels = newChoice; // } @@ -182,6 +182,46 @@ namespace storm { // return storm::models::Dtmc<ValueType>(newMatBuilder.build(), newLabeling, newStateRewards, std::move(newTransitionRewards), newChoiceLabels); } + template<typename ValueType> + Dtmc<ValueType>::ConstraintCollector::ConstraintCollector(Dtmc<ValueType> const& dtmc) { + process(dtmc); + } + + template<typename ValueType> + std::unordered_set<carl::Constraint<ValueType>> const& Dtmc<ValueType>::ConstraintCollector::getWellformedConstraints() const { + return this->wellformedConstraintSet; + } + + template<typename ValueType> + std::unordered_set<carl::Constraint<ValueType>> const& Dtmc<ValueType>::ConstraintCollector::getGraphPreservingConstraints() const { + return this->graphPreservingConstraintSet; + } + + template<typename ValueType> + void Dtmc<ValueType>::ConstraintCollector::process(storm::models::sparse::Dtmc<ValueType> const& dtmc) { + for(uint_fast64_t state = 0; state < dtmc.getNumberOfStates(); ++state) { + ValueType sum = storm::utility::zero<ValueType>(); + for (auto const& transition : dtmc.getRows(state)) { + sum += transition.getValue(); + if (!comparator.isConstant(transition.getValue())) { + wellformedConstraintSet.emplace(transition.getValue() - 1, storm::CompareRelation::LEQ); + wellformedConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GEQ); + graphPreservingConstraintSet.emplace(transition.getValue(), storm::CompareRelation::GT); + } + } + STORM_LOG_ASSERT(!comparator.isConstant(sum) || comparator.isOne(sum), "If the sum is a constant, it must be equal to 1."); + if(!comparator.isConstant(sum)) { + wellformedConstraintSet.emplace(sum - 1, storm::CompareRelation::EQ); + } + + } + } + + template<typename ValueType> + void Dtmc<ValueType>::ConstraintCollector::operator()(storm::models::sparse::Dtmc<ValueType> const& dtmc) { + process(dtmc); + } + template <typename ValueType> bool Dtmc<ValueType>::checkValidityOfProbabilityMatrix() const { if (this->getTransitionMatrix().getRowCount() != this->getTransitionMatrix().getColumnCount()) { diff --git a/src/models/sparse/Dtmc.h b/src/models/sparse/Dtmc.h index 2ff6dd884..37095f05e 100644 --- a/src/models/sparse/Dtmc.h +++ b/src/models/sparse/Dtmc.h @@ -27,7 +27,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -41,7 +41,7 @@ namespace storm { Dtmc(storm::storage::SparseMatrix<ValueType>&& probabilityMatrix, storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); Dtmc(Dtmc<ValueType> const& dtmc) = default; Dtmc& operator=(Dtmc<ValueType> const& dtmc) = default; @@ -59,6 +59,57 @@ namespace storm { */ Dtmc<ValueType> getSubDtmc(storm::storage::BitVector const& states) const; + class ConstraintCollector { + private: + // A set of constraints that says that the DTMC actually has valid probability distributions in all states. + std::unordered_set<carl::Constraint<ValueType>> wellformedConstraintSet; + + // A set of constraints that makes sure that the underlying graph of the model does not change depending + // on the parameter values. + std::unordered_set<carl::Constraint<ValueType>> graphPreservingConstraintSet; + + // A comparator that is used for + storm::utility::ConstantsComparator<ValueType> comparator; + + public: + /*! + * Constructs the a constraint collector for the given DTMC. The constraints are built and ready for + * retrieval after the construction. + * + * @param dtmc The DTMC for which to create the constraints. + */ + ConstraintCollector(storm::models::sparse::Dtmc<ValueType> const& dtmc); + + /*! + * Returns the set of wellformed-ness constraints. + * + * @return The set of wellformed-ness constraints. + */ + std::unordered_set<carl::Constraint<ValueType>> const& getWellformedConstraints() const; + + /*! + * Returns the set of graph-preserving constraints. + * + * @return The set of graph-preserving constraints. + */ + std::unordered_set<carl::Constraint<ValueType>> const& getGraphPreservingConstraints() const; + + /*! + * Constructs the constraints for the given DTMC. + * + * @param dtmc The DTMC for which to create the constraints. + */ + void process(storm::models::sparse::Dtmc<ValueType> const& dtmc); + + /*! + * Constructs the constraints for the given DTMC by calling the process method. + * + * @param dtmc The DTMC for which to create the constraints. + */ + void operator()(storm::models::sparse::Dtmc<ValueType> const& dtmc); + + }; + private: /*! * Checks the probability matrix for validity. diff --git a/src/models/sparse/MarkovAutomaton.cpp b/src/models/sparse/MarkovAutomaton.cpp index 4032452d8..d161a64f6 100644 --- a/src/models/sparse/MarkovAutomaton.cpp +++ b/src/models/sparse/MarkovAutomaton.cpp @@ -13,7 +13,7 @@ namespace storm { std::vector<ValueType> const& exitRates, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : NondeterministicModel<ValueType>(storm::models::ModelType::MarkovAutomaton, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling), markovianStates(markovianStates), exitRates(exitRates), closed(false) { this->turnRatesToProbabilities(); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "There are transition rewards for nonexistent transitions."); @@ -26,7 +26,7 @@ namespace storm { std::vector<ValueType> const& exitRates, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : NondeterministicModel<ValueType>(storm::models::ModelType::MarkovAutomaton, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)), markovianStates(markovianStates), exitRates(std::move(exitRates)), closed(false) { this->turnRatesToProbabilities(); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "There are transition rewards for nonexistent transitions."); diff --git a/src/models/sparse/MarkovAutomaton.h b/src/models/sparse/MarkovAutomaton.h index d051965cc..622da37b0 100644 --- a/src/models/sparse/MarkovAutomaton.h +++ b/src/models/sparse/MarkovAutomaton.h @@ -31,7 +31,7 @@ namespace storm { std::vector<ValueType> const& exitRates, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -50,7 +50,7 @@ namespace storm { std::vector<ValueType> const& exitRates, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); MarkovAutomaton(MarkovAutomaton const& other) = default; MarkovAutomaton& operator=(MarkovAutomaton const& other) = default; diff --git a/src/models/sparse/Mdp.cpp b/src/models/sparse/Mdp.cpp index f8775107b..1c8bd1b50 100644 --- a/src/models/sparse/Mdp.cpp +++ b/src/models/sparse/Mdp.cpp @@ -11,7 +11,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : NondeterministicModel<ValueType>(storm::models::ModelType::Mdp, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid."); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist."); @@ -23,20 +23,20 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : NondeterministicModel<ValueType>(storm::models::ModelType::Mdp, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { STORM_LOG_THROW(this->checkValidityOfProbabilityMatrix(), storm::exceptions::InvalidArgumentException, "The probability matrix is invalid."); STORM_LOG_THROW(!this->hasTransitionRewards() || this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), storm::exceptions::InvalidArgumentException, "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist."); } template <typename ValueType> - Mdp<ValueType> Mdp<ValueType>::restrictChoiceLabels(boost::container::flat_set<uint_fast64_t> const& enabledChoiceLabels) const { + Mdp<ValueType> Mdp<ValueType>::restrictChoiceLabels(LabelSet const& enabledChoiceLabels) const { STORM_LOG_THROW(this->hasChoiceLabeling(), storm::exceptions::InvalidArgumentException, "Restriction to label set is impossible for unlabeled model."); - std::vector<boost::container::flat_set<uint_fast64_t>> const& choiceLabeling = this->getChoiceLabeling(); + std::vector<LabelSet> const& choiceLabeling = this->getChoiceLabeling(); storm::storage::SparseMatrixBuilder<ValueType> transitionMatrixBuilder(0, this->getTransitionMatrix().getColumnCount(), 0, true, true); - std::vector<boost::container::flat_set<uint_fast64_t>> newChoiceLabeling; + std::vector<LabelSet> newChoiceLabeling; // Check for each choice of each state, whether the choice labels are fully contained in the given label set. uint_fast64_t currentRow = 0; @@ -71,7 +71,7 @@ namespace storm { Mdp<ValueType> restrictedMdp(transitionMatrixBuilder.build(), storm::models::sparse::StateLabeling(this->getStateLabeling()), this->hasStateRewards() ? boost::optional<std::vector<ValueType>>(this->getStateRewardVector()) : boost::optional<std::vector<ValueType>>(), this->hasTransitionRewards() ? boost::optional<storm::storage::SparseMatrix<ValueType>>(this->getTransitionRewardMatrix()) : boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>(newChoiceLabeling)); + boost::optional<std::vector<LabelSet>>(newChoiceLabeling)); return restrictedMdp; } diff --git a/src/models/sparse/Mdp.h b/src/models/sparse/Mdp.h index 4b62f7377..18879ee76 100644 --- a/src/models/sparse/Mdp.h +++ b/src/models/sparse/Mdp.h @@ -27,7 +27,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -42,7 +42,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); Mdp(Mdp const& other) = default; Mdp& operator=(Mdp const& other) = default; @@ -61,7 +61,7 @@ namespace storm { * and which ones need to be ignored. * @return A restricted version of the current MDP that only uses choice labels from the given set. */ - Mdp<ValueType> restrictChoiceLabels(boost::container::flat_set<uint_fast64_t> const& enabledChoiceLabels) const; + Mdp<ValueType> restrictChoiceLabels(LabelSet const& enabledChoiceLabels) const; private: /*! diff --git a/src/models/sparse/Model.cpp b/src/models/sparse/Model.cpp index e1230d051..b2ca03c1d 100644 --- a/src/models/sparse/Model.cpp +++ b/src/models/sparse/Model.cpp @@ -12,7 +12,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : ModelBase(modelType), transitionMatrix(transitionMatrix), stateLabeling(stateLabeling), @@ -28,7 +28,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : ModelBase(modelType), transitionMatrix(std::move(transitionMatrix)), stateLabeling(std::move(stateLabeling)), @@ -109,7 +109,7 @@ namespace storm { } template <typename ValueType> - std::vector<boost::container::flat_set<uint_fast64_t>> const& Model<ValueType>::getChoiceLabeling() const { + std::vector<LabelSet> const& Model<ValueType>::getChoiceLabeling() const { return choiceLabeling.get(); } @@ -159,7 +159,7 @@ namespace storm { result += getTransitionRewardMatrix().getSizeInBytes(); } if (hasChoiceLabeling()) { - result += getChoiceLabeling().size() * sizeof(boost::container::flat_set<uint_fast64_t>); + result += getChoiceLabeling().size() * sizeof(LabelSet); } return result; } diff --git a/src/models/sparse/Model.h b/src/models/sparse/Model.h index 103da9fba..5457370e6 100644 --- a/src/models/sparse/Model.h +++ b/src/models/sparse/Model.h @@ -16,6 +16,9 @@ namespace storm { namespace models { namespace sparse { + // The type used for storing a set of labels. + typedef boost::container::flat_set<uint_fast64_t> LabelSet; + /*! * Base class for all sparse models. */ @@ -45,7 +48,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -62,7 +65,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Retrieves the backward transition relation of the model, i.e. a set of transitions between states @@ -175,7 +178,7 @@ namespace storm { * * @return The labels for the choices of the model. */ - std::vector<boost::container::flat_set<uint_fast64_t>> const& getChoiceLabeling() const; + std::vector<LabelSet> const& getChoiceLabeling() const; /*! * Returns the state labeling associated with this model. @@ -287,7 +290,7 @@ namespace storm { boost::optional<storm::storage::SparseMatrix<ValueType>> transitionRewardMatrix; // If set, a vector representing the labels of choices. - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> choiceLabeling; + boost::optional<std::vector<LabelSet>> choiceLabeling; }; } // namespace sparse diff --git a/src/models/sparse/NondeterministicModel.cpp b/src/models/sparse/NondeterministicModel.cpp index 860418f1e..57dbe4b84 100644 --- a/src/models/sparse/NondeterministicModel.cpp +++ b/src/models/sparse/NondeterministicModel.cpp @@ -12,7 +12,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling) : Model<ValueType>(modelType, transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { // Intentionally left empty. } @@ -23,7 +23,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix, - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling) + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling) : Model<ValueType>(modelType, std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { // Intentionally left empty. diff --git a/src/models/sparse/NondeterministicModel.h b/src/models/sparse/NondeterministicModel.h index 9151c2145..27cdaf11e 100644 --- a/src/models/sparse/NondeterministicModel.h +++ b/src/models/sparse/NondeterministicModel.h @@ -29,7 +29,7 @@ namespace storm { storm::models::sparse::StateLabeling const& stateLabeling, boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>> const& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> const& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); /*! * Constructs a model by moving the given data. @@ -46,7 +46,7 @@ namespace storm { storm::models::sparse::StateLabeling&& stateLabeling, boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), boost::optional<storm::storage::SparseMatrix<ValueType>>&& optionalTransitionRewardMatrix = boost::optional<storm::storage::SparseMatrix<ValueType>>(), - boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>&& optionalChoiceLabeling = boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>()); + boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>()); NondeterministicModel(NondeterministicModel const& other) = default; NondeterministicModel& operator=(NondeterministicModel const& other) = default; diff --git a/src/models/sparse/StochasticTwoPlayerGame.cpp b/src/models/sparse/StochasticTwoPlayerGame.cpp new file mode 100644 index 000000000..8e6daa52a --- /dev/null +++ b/src/models/sparse/StochasticTwoPlayerGame.cpp @@ -0,0 +1,41 @@ +#include "src/models/sparse/StochasticTwoPlayerGame.h" + +#include "src/adapters/CarlAdapter.h" + +namespace storm { + namespace models { + namespace sparse { + + template <typename ValueType> + StochasticTwoPlayerGame<ValueType>::StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix, + storm::storage::SparseMatrix<ValueType> const& player2Matrix, + storm::models::sparse::StateLabeling const& stateLabeling, + boost::optional<std::vector<ValueType>> const& optionalStateRewardVector, + boost::optional<std::vector<LabelSet>> const& optionalPlayer1ChoiceLabeling, + boost::optional<std::vector<LabelSet>> const& optionalPlayer2ChoiceLabeling) + : NondeterministicModel<ValueType>(storm::models::ModelType::S2pg, player2Matrix, stateLabeling, optionalStateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>>(), optionalPlayer2ChoiceLabeling), player1Matrix(player1Matrix), player1Labels(optionalPlayer1ChoiceLabeling) { + // Intentionally left empty. + } + + + template <typename ValueType> + StochasticTwoPlayerGame<ValueType>::StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type>&& player1Matrix, + storm::storage::SparseMatrix<ValueType>&& player2Matrix, + storm::models::sparse::StateLabeling&& stateLabeling, + boost::optional<std::vector<ValueType>>&& optionalStateRewardVector, + boost::optional<std::vector<LabelSet>>&& optionalPlayer1ChoiceLabeling, + boost::optional<std::vector<LabelSet>>&& optionalPlayer2ChoiceLabeling) + : NondeterministicModel<ValueType>(storm::models::ModelType::S2pg, std::move(player2Matrix), std::move(stateLabeling), std::move(optionalStateRewardVector), boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(optionalPlayer2ChoiceLabeling)), player1Matrix(std::move(player1Matrix)), player1Labels(std::move(optionalPlayer1ChoiceLabeling)) { + // Intentionally left empty. + } + + template class StochasticTwoPlayerGame<double>; + template class StochasticTwoPlayerGame<float>; + +#ifdef STORM_HAVE_CARL + template class StochasticTwoPlayerGame<storm::RationalFunction>; +#endif + + } // namespace sparse + } // namespace models +} // namespace storm \ No newline at end of file diff --git a/src/models/sparse/StochasticTwoPlayerGame.h b/src/models/sparse/StochasticTwoPlayerGame.h new file mode 100644 index 000000000..cbd2988ca --- /dev/null +++ b/src/models/sparse/StochasticTwoPlayerGame.h @@ -0,0 +1,77 @@ +#ifndef STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_ +#define STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_ + +#include "src/models/sparse/NondeterministicModel.h" +#include "src/utility/OsDetection.h" + +namespace storm { + namespace models { + namespace sparse { + + /*! + * This class represents a (discrete-time) stochastic two-player game. + */ + template <typename ValueType> + class StochasticTwoPlayerGame : public NondeterministicModel<ValueType> { + public: + /*! + * Constructs a model from the given data. + * + * @param player1Matrix The matrix representing the choices of player 1. + * @param player2Matrix The matrix representing the choices of player 2. + * @param stateLabeling The labeling of the states. + * @param optionalStateRewardVector The reward values associated with the states. + * @param optionalPlayer1ChoiceLabeling A vector that represents the labels associated with the choices of each player 1 state. + * @param optionalPlayer2ChoiceLabeling A vector that represents the labels associated with the choices of each player 2 state. + */ + StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix, + storm::storage::SparseMatrix<ValueType> const& player2Matrix, + storm::models::sparse::StateLabeling const& stateLabeling, + boost::optional<std::vector<ValueType>> const& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), + boost::optional<std::vector<LabelSet>> const& optionalPlayer1ChoiceLabeling = boost::optional<std::vector<LabelSet>>(), + boost::optional<std::vector<LabelSet>> const& optionalPlayer2ChoiceLabeling = boost::optional<std::vector<LabelSet>>()); + + /*! + * Constructs a model by moving the given data. + * + * @param player1Matrix The matrix representing the choices of player 1. + * @param player2Matrix The matrix representing the choices of player 2. + * @param stateLabeling The labeling of the states. + * @param optionalStateRewardVector The reward values associated with the states. + * @param optionalPlayer1ChoiceLabeling A vector that represents the labels associated with the choices of each player 1 state. + * @param optionalPlayer2ChoiceLabeling A vector that represents the labels associated with the choices of each player 2 state. + */ + StochasticTwoPlayerGame(storm::storage::SparseMatrix<storm::storage::sparse::state_type>&& player1Matrix, + storm::storage::SparseMatrix<ValueType>&& player2Matrix, + storm::models::sparse::StateLabeling&& stateLabeling, + boost::optional<std::vector<ValueType>>&& optionalStateRewardVector = boost::optional<std::vector<ValueType>>(), + boost::optional<std::vector<LabelSet>>&& optionalPlayer1ChoiceLabeling = boost::optional<std::vector<LabelSet>>(), + boost::optional<std::vector<LabelSet>>&& optionalPlayer2ChoiceLabeling = boost::optional<std::vector<LabelSet>>()); + + StochasticTwoPlayerGame(StochasticTwoPlayerGame const& other) = default; + StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame const& other) = default; + +#ifndef WINDOWS + StochasticTwoPlayerGame(StochasticTwoPlayerGame&& other) = default; + StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame&& other) = default; +#endif + + private: + // A matrix that stores the player 1 choices. This matrix contains a row group for each player 1 node. Every + // row group contains a row for each choice in that player 1 node. Each such row contains exactly one + // (non-zero) entry at a column that indicates the player 2 node this choice leads to (which is essentially + // the index of a row group in the matrix for player 2). + storm::storage::SparseMatrix<storm::storage::sparse::state_type> player1Matrix; + + // An (optional) vector of labels attached to the choices of player 1. Each row of the matrix can be equipped + // with a set of labels to tag certain choices. + boost::optional<std::vector<LabelSet>> player1Labels; + + // The matrix and labels for player 2 are stored in the superclass. + }; + + } // namespace sparse + } // namespace models +} // namespace storm + +#endif /* STORM_MODELS_SPARSE_STOCHASTICTWOPLAYERGAME_H_ */ diff --git a/src/models/symbolic/StochasticTwoPlayerGame.cpp b/src/models/symbolic/StochasticTwoPlayerGame.cpp new file mode 100644 index 000000000..0424f100c --- /dev/null +++ b/src/models/symbolic/StochasticTwoPlayerGame.cpp @@ -0,0 +1,42 @@ +#include "src/models/symbolic/StochasticTwoPlayerGame.h" + +namespace storm { + namespace models { + namespace symbolic { + + template<storm::dd::DdType Type> + StochasticTwoPlayerGame<Type>::StochasticTwoPlayerGame(std::shared_ptr<storm::dd::DdManager<Type>> manager, + storm::dd::Bdd<Type> reachableStates, + storm::dd::Bdd<Type> initialStates, + storm::dd::Add<Type> transitionMatrix, + std::set<storm::expressions::Variable> const& rowVariables, + std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> rowExpressionAdapter, + std::set<storm::expressions::Variable> const& columnVariables, + std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> columnExpressionAdapter, + std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, + std::set<storm::expressions::Variable> const& player1Variables, + std::set<storm::expressions::Variable> const& player2Variables, + std::set<storm::expressions::Variable> const& nondeterminismVariables, + std::map<std::string, storm::expressions::Expression> labelToExpressionMap, + boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector, + boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix) + : NondeterministicModel<Type>(storm::models::ModelType::S2pg, manager, reachableStates, initialStates, transitionMatrix, rowVariables, rowExpressionAdapter, columnVariables, columnExpressionAdapter, rowColumnMetaVariablePairs, nondeterminismVariables, labelToExpressionMap, optionalStateRewardVector, optionalTransitionRewardMatrix), player1Variables(player1Variables), player2Variables(player2Variables) { + // Intentionally left empty. + } + + template<storm::dd::DdType Type> + std::set<storm::expressions::Variable> const& StochasticTwoPlayerGame<Type>::getPlayer1Variables() const { + return player1Variables; + } + + template<storm::dd::DdType Type> + std::set<storm::expressions::Variable> const& StochasticTwoPlayerGame<Type>::getPlayer2Variables() const { + return player2Variables; + } + + // Explicitly instantiate the template class. + template class StochasticTwoPlayerGame<storm::dd::DdType::CUDD>; + + } // namespace symbolic + } // namespace models +} // namespace storm \ No newline at end of file diff --git a/src/models/symbolic/StochasticTwoPlayerGame.h b/src/models/symbolic/StochasticTwoPlayerGame.h new file mode 100644 index 000000000..841cce5c4 --- /dev/null +++ b/src/models/symbolic/StochasticTwoPlayerGame.h @@ -0,0 +1,88 @@ +#ifndef STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_ +#define STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_ + +#include "src/models/symbolic/NondeterministicModel.h" +#include "src/utility/OsDetection.h" + +namespace storm { + namespace models { + namespace symbolic { + + /*! + * This class represents a discrete-time stochastic two-player game. + */ + template<storm::dd::DdType Type> + class StochasticTwoPlayerGame : public NondeterministicModel<Type> { + public: + StochasticTwoPlayerGame(StochasticTwoPlayerGame<Type> const& other) = default; + StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame<Type> const& other) = default; + +#ifndef WINDOWS + StochasticTwoPlayerGame(StochasticTwoPlayerGame<Type>&& other) = default; + StochasticTwoPlayerGame& operator=(StochasticTwoPlayerGame<Type>&& other) = default; +#endif + + /*! + * Constructs a model from the given data. + * + * @param manager The manager responsible for the decision diagrams. + * @param reachableStates A DD representing the reachable states. + * @param initialStates A DD representing the initial states of the model. + * @param transitionMatrix The matrix representing the transitions in the model. + * @param rowVariables The set of row meta variables used in the DDs. + * @param rowExpressionAdapter An object that can be used to translate expressions in terms of the row + * meta variables. + * @param columVariables The set of column meta variables used in the DDs. + * @param columnExpressionAdapter An object that can be used to translate expressions in terms of the + * column meta variables. + * @param rowColumnMetaVariablePairs All pairs of row/column meta variables. + * @param player1Variables The meta variables used to encode the nondeterministic choices of player 1. + * @param player2Variables The meta variables used to encode the nondeterministic choices of player 2. + * @param allNondeterminismVariables The meta variables used to encode the nondeterminism in the model. + * @param labelToExpressionMap A mapping from label names to their defining expressions. + * @param optionalStateRewardVector The reward values associated with the states. + * @param optionalTransitionRewardMatrix The reward values associated with the transitions of the model. + */ + StochasticTwoPlayerGame(std::shared_ptr<storm::dd::DdManager<Type>> manager, + storm::dd::Bdd<Type> reachableStates, + storm::dd::Bdd<Type> initialStates, + storm::dd::Add<Type> transitionMatrix, + std::set<storm::expressions::Variable> const& rowVariables, + std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> rowExpressionAdapter, + std::set<storm::expressions::Variable> const& columnVariables, + std::shared_ptr<storm::adapters::AddExpressionAdapter<Type>> columnExpressionAdapter, + std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, + std::set<storm::expressions::Variable> const& player1Variables, + std::set<storm::expressions::Variable> const& player2Variables, + std::set<storm::expressions::Variable> const& allNondeterminismVariables, + std::map<std::string, storm::expressions::Expression> labelToExpressionMap = std::map<std::string, storm::expressions::Expression>(), + boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector = boost::optional<storm::dd::Dd<Type>>(), + boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix = boost::optional<storm::dd::Dd<Type>>()); + + /*! + * Retrieeves the set of meta variables used to encode the nondeterministic choices of player 1. + * + * @return The set of meta variables used to encode the nondeterministic choices of player 1. + */ + std::set<storm::expressions::Variable> const& getPlayer1Variables() const; + + /*! + * Retrieeves the set of meta variables used to encode the nondeterministic choices of player 2. + * + * @return The set of meta variables used to encode the nondeterministic choices of player 2. + */ + std::set<storm::expressions::Variable> const& getPlayer2Variables() const; + + private: + // The meta variables used to encode the nondeterministic choices of player 1. + std::set<storm::expressions::Variable> player1Variables; + + // The meta variables used to encode the nondeterministic choices of player 2. + std::set<storm::expressions::Variable> player2Variables; + }; + + } // namespace symbolic + } // namespace models +} // namespace storm + +#endif /* STORM_MODELS_SYMBOLIC_STOCHASTICTWOPLAYERGAME_H_ */ diff --git a/src/parser/ExpressionParser.cpp b/src/parser/ExpressionParser.cpp index 5666b39db..507753a71 100644 --- a/src/parser/ExpressionParser.cpp +++ b/src/parser/ExpressionParser.cpp @@ -23,7 +23,7 @@ namespace storm { } minMaxExpression.name("min/max expression"); - identifierExpression = identifier[qi::_val = phoenix::bind(&ExpressionParser::getIdentifierExpression, phoenix::ref(*this), qi::_1)]; + identifierExpression = identifier[qi::_val = phoenix::bind(&ExpressionParser::getIdentifierExpression, phoenix::ref(*this), qi::_1, allowBacktracking, phoenix::ref(qi::_pass))]; identifierExpression.name("identifier expression"); literalExpression = trueFalse_[qi::_val = qi::_1] | strict_double[qi::_val = phoenix::bind(&ExpressionParser::createDoubleLiteralExpression, phoenix::ref(*this), qi::_1, qi::_pass)] | qi::int_[qi::_val = phoenix::bind(&ExpressionParser::createIntegerLiteralExpression, phoenix::ref(*this), qi::_1)]; @@ -53,7 +53,7 @@ namespace storm { plusExpression.name("plus expression"); if (allowBacktracking) { - relativeExpression = plusExpression[qi::_val = qi::_1] > -(relationalOperator_ >> plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; + relativeExpression = plusExpression[qi::_val = qi::_1] >> -(relationalOperator_ >> plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; } else { relativeExpression = plusExpression[qi::_val = qi::_1] > -(relationalOperator_ > plusExpression)[qi::_val = phoenix::bind(&ExpressionParser::createRelationalExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; } @@ -70,7 +70,7 @@ namespace storm { andExpression.name("and expression"); if (allowBacktracking) { - orExpression = andExpression[qi::_val = qi::_1] > *(orOperator_ >> andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; + orExpression = andExpression[qi::_val = qi::_1] >> *(orOperator_ >> andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; } else { orExpression = andExpression[qi::_val = qi::_1] > *(orOperator_ > andExpression)[qi::_val = phoenix::bind(&ExpressionParser::createOrExpression, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)]; } @@ -82,6 +82,24 @@ namespace storm { expression %= iteExpression; expression.name("expression"); + /*! + * Enable this for debugging purposes. + debug(expression); + debug(iteExpression); + debug(orExpression); + debug(andExpression); + debug(equalityExpression); + debug(relativeExpression); + debug(plusExpression); + debug(multiplicationExpression); + debug(powerExpression); + debug(unaryExpression); + debug(atomicExpression); + debug(literalExpression); + debug(identifierExpression); + */ + + // Enable error reporting. qi::on_error<qi::fail>(expression, handler(qi::_1, qi::_2, qi::_3, qi::_4)); qi::on_error<qi::fail>(iteExpression, handler(qi::_1, qi::_2, qi::_3, qi::_4)); @@ -304,11 +322,18 @@ namespace storm { return manager.boolean(false); } - storm::expressions::Expression ExpressionParser::getIdentifierExpression(std::string const& identifier) const { + storm::expressions::Expression ExpressionParser::getIdentifierExpression(std::string const& identifier, bool allowBacktracking, bool& pass) const { if (this->createExpressions) { STORM_LOG_THROW(this->identifiers_ != nullptr, storm::exceptions::WrongFormatException, "Unable to substitute identifier expressions without given mapping."); storm::expressions::Expression const* expression = this->identifiers_->find(identifier); - STORM_LOG_THROW(expression != nullptr, storm::exceptions::WrongFormatException, "Parsing error in line " << get_line(qi::_3) << ": Undeclared identifier '" << identifier << "'."); + if (expression == nullptr) { + if (allowBacktracking) { + pass = false; + return manager.boolean(false); + } else { + STORM_LOG_THROW(expression != nullptr, storm::exceptions::WrongFormatException, "Parsing error in line " << get_line(qi::_3) << ": Undeclared identifier '" << identifier << "'."); + } + } return *expression; } else { return manager.boolean(false); diff --git a/src/parser/ExpressionParser.h b/src/parser/ExpressionParser.h index 2a67e7df6..55df73e2a 100644 --- a/src/parser/ExpressionParser.h +++ b/src/parser/ExpressionParser.h @@ -223,7 +223,7 @@ namespace storm { storm::expressions::Expression createIntegerLiteralExpression(int value) const; storm::expressions::Expression createMinimumMaximumExpression(storm::expressions::Expression const& e1, storm::expressions::OperatorType const& operatorType, storm::expressions::Expression const& e2) const; storm::expressions::Expression createFloorCeilExpression(storm::expressions::OperatorType const& operatorType, storm::expressions::Expression const& e1) const; - storm::expressions::Expression getIdentifierExpression(std::string const& identifier) const; + storm::expressions::Expression getIdentifierExpression(std::string const& identifier, bool allowBacktracking, bool& pass) const; bool isValidIdentifier(std::string const& identifier); diff --git a/src/parser/FormulaParser.cpp b/src/parser/FormulaParser.cpp index 218a51734..268865f47 100644 --- a/src/parser/FormulaParser.cpp +++ b/src/parser/FormulaParser.cpp @@ -38,7 +38,10 @@ namespace storm { booleanLiteralFormula = (qi::lit("true")[qi::_a = true] | qi::lit("false")[qi::_a = false])[qi::_val = phoenix::bind(&FormulaParser::createBooleanLiteralFormula, phoenix::ref(*this), qi::_a)]; booleanLiteralFormula.name("boolean literal formula"); - atomicStateFormula = booleanLiteralFormula | labelFormula | expressionFormula | (qi::lit("(") > stateFormula > qi::lit(")")); + operatorFormula = probabilityOperator | rewardOperator | steadyStateOperator; + operatorFormula.name("operator formulas"); + + atomicStateFormula = booleanLiteralFormula | labelFormula | expressionFormula | (qi::lit("(") > stateFormula > qi::lit(")")) | operatorFormula; atomicStateFormula.name("atomic state formula"); notStateFormula = (-unaryBooleanOperator_ >> atomicStateFormula)[qi::_val = phoenix::bind(&FormulaParser::createUnaryBooleanStateFormula, phoenix::ref(*this), qi::_2, qi::_1)]; @@ -89,7 +92,7 @@ namespace storm { orStateFormula = andStateFormula[qi::_val = qi::_1] >> *(qi::lit("|") >> andStateFormula)[qi::_val = phoenix::bind(&FormulaParser::createBinaryBooleanStateFormula, phoenix::ref(*this), qi::_val, qi::_1, storm::logic::BinaryBooleanStateFormula::OperatorType::Or)]; orStateFormula.name("or state formula"); - stateFormula = (probabilityOperator | rewardOperator | steadyStateOperator | orStateFormula); + stateFormula = (orStateFormula); stateFormula.name("state formula"); start = qi::eps > stateFormula >> qi::skip(boost::spirit::ascii::space | qi::lit("//") >> *(qi::char_ - (qi::eol | qi::eoi)))[qi::eps] >> qi::eoi; diff --git a/src/parser/FormulaParser.h b/src/parser/FormulaParser.h index c5735286b..d2de1a73f 100644 --- a/src/parser/FormulaParser.h +++ b/src/parser/FormulaParser.h @@ -132,6 +132,7 @@ namespace storm { qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> pathFormulaWithoutUntil; qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> simplePathFormula; qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> atomicStateFormula; + qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> operatorFormula; qi::rule<Iterator, std::string(), Skipper> label; qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> andStateFormula; diff --git a/src/settings/modules/GmmxxEquationSolverSettings.cpp b/src/settings/modules/GmmxxEquationSolverSettings.cpp index 4b9f27f9e..359836968 100644 --- a/src/settings/modules/GmmxxEquationSolverSettings.cpp +++ b/src/settings/modules/GmmxxEquationSolverSettings.cpp @@ -7,7 +7,7 @@ namespace storm { namespace modules { const std::string GmmxxEquationSolverSettings::moduleName = "gmm++"; - const std::string GmmxxEquationSolverSettings::techniqueOptionName = "tech"; + const std::string GmmxxEquationSolverSettings::techniqueOptionName = "method"; const std::string GmmxxEquationSolverSettings::preconditionOptionName = "precond"; const std::string GmmxxEquationSolverSettings::restartOptionName = "restart"; const std::string GmmxxEquationSolverSettings::maximalIterationsOptionName = "maxiter"; @@ -32,38 +32,38 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); } - bool GmmxxEquationSolverSettings::isLinearEquationSystemTechniqueSet() const { + bool GmmxxEquationSolverSettings::isLinearEquationSystemMethodSet() const { return this->getOption(techniqueOptionName).getHasOptionBeenSet(); } - GmmxxEquationSolverSettings::LinearEquationTechnique GmmxxEquationSolverSettings::getLinearEquationSystemTechnique() const { + GmmxxEquationSolverSettings::LinearEquationMethod GmmxxEquationSolverSettings::getLinearEquationSystemMethod() const { std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString(); if (linearEquationSystemTechniqueAsString == "bicgstab") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab; + return GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab; } else if (linearEquationSystemTechniqueAsString == "qmr") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr; + return GmmxxEquationSolverSettings::LinearEquationMethod::Qmr; } else if (linearEquationSystemTechniqueAsString == "gmres") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres; + return GmmxxEquationSolverSettings::LinearEquationMethod::Gmres; } else if (linearEquationSystemTechniqueAsString == "jacobi") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi; + return GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } - bool GmmxxEquationSolverSettings::isPreconditioningTechniqueSet() const { + bool GmmxxEquationSolverSettings::isPreconditioningMethodSet() const { return this->getOption(preconditionOptionName).getHasOptionBeenSet(); } - GmmxxEquationSolverSettings::PreconditioningTechnique GmmxxEquationSolverSettings::getPreconditioningTechnique() const { - std::string preconditioningTechniqueAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString(); - if (preconditioningTechniqueAsString == "ilu") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu; - } else if (preconditioningTechniqueAsString == "diagonal") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal; - } else if (preconditioningTechniqueAsString == "none") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::None; + GmmxxEquationSolverSettings::PreconditioningMethod GmmxxEquationSolverSettings::getPreconditioningMethod() const { + std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString(); + if (PreconditioningMethodAsString == "ilu") { + return GmmxxEquationSolverSettings::PreconditioningMethod::Ilu; + } else if (PreconditioningMethodAsString == "diagonal") { + return GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal; + } else if (PreconditioningMethodAsString == "none") { + return GmmxxEquationSolverSettings::PreconditioningMethod::None; } - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << preconditioningTechniqueAsString << "' selected."); + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected."); } bool GmmxxEquationSolverSettings::isRestartIterationCountSet() const { @@ -99,7 +99,7 @@ namespace storm { } bool GmmxxEquationSolverSettings::check() const { - bool optionsSet = isLinearEquationSystemTechniqueSet() || isPreconditioningTechniqueSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet(); + bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet(); STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect."); diff --git a/src/settings/modules/GmmxxEquationSolverSettings.h b/src/settings/modules/GmmxxEquationSolverSettings.h index 07bc90174..8cb0761a8 100644 --- a/src/settings/modules/GmmxxEquationSolverSettings.h +++ b/src/settings/modules/GmmxxEquationSolverSettings.h @@ -12,11 +12,11 @@ namespace storm { */ class GmmxxEquationSolverSettings : public ModuleSettings { public: - // An enumeration of all available techniques for solving linear equations. - enum class LinearEquationTechnique { Bicgstab, Qmr, Gmres, Jacobi }; + // An enumeration of all available methods for solving linear equations. + enum class LinearEquationMethod { Bicgstab, Qmr, Gmres, Jacobi }; - // An enumeration of all available preconditioning techniques. - enum class PreconditioningTechnique { Ilu, Diagonal, None }; + // An enumeration of all available preconditioning methods. + enum class PreconditioningMethod { Ilu, Diagonal, None }; // An enumeration of all available convergence criteria. enum class ConvergenceCriterion { Absolute, Relative }; @@ -29,32 +29,32 @@ namespace storm { GmmxxEquationSolverSettings(storm::settings::SettingsManager& settingsManager); /*! - * Retrieves whether the linear equation system technique has been set. + * Retrieves whether the linear equation system method has been set. * - * @return True iff the linear equation system technique has been set. + * @return True iff the linear equation system method has been set. */ - bool isLinearEquationSystemTechniqueSet() const; + bool isLinearEquationSystemMethodSet() const; /*! - * Retrieves the technique that is to be used for solving systems of linear equations. + * Retrieves the method that is to be used for solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - LinearEquationTechnique getLinearEquationSystemTechnique() const; + LinearEquationMethod getLinearEquationSystemMethod() const; /*! - * Retrieves whether the preconditioning technique has been set. + * Retrieves whether the preconditioning method has been set. * - * @return True iff the preconditioning technique has been set. + * @return True iff the preconditioning method has been set. */ - bool isPreconditioningTechniqueSet() const; + bool isPreconditioningMethodSet() const; /*! - * Retrieves the technique that is to be used for preconditioning solving systems of linear equations. + * Retrieves the method that is to be used for preconditioning solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - PreconditioningTechnique getPreconditioningTechnique() const; + PreconditioningMethod getPreconditioningMethod() const; /*! * Retrieves whether the restart iteration count has been set. @@ -64,7 +64,7 @@ namespace storm { bool isRestartIterationCountSet() const; /*! - * Retrieves the number of iterations after which restarted techniques are to be restarted. + * Retrieves the number of iterations after which restarted methods are to be restarted. * * @return The number of iterations after which to restart. */ diff --git a/src/settings/modules/NativeEquationSolverSettings.cpp b/src/settings/modules/NativeEquationSolverSettings.cpp index 76af99b3e..482713ed8 100644 --- a/src/settings/modules/NativeEquationSolverSettings.cpp +++ b/src/settings/modules/NativeEquationSolverSettings.cpp @@ -7,20 +7,23 @@ namespace storm { namespace modules { const std::string NativeEquationSolverSettings::moduleName = "native"; - const std::string NativeEquationSolverSettings::techniqueOptionName = "tech"; + const std::string NativeEquationSolverSettings::techniqueOptionName = "method"; + const std::string NativeEquationSolverSettings::omegaOptionName = "soromega"; const std::string NativeEquationSolverSettings::maximalIterationsOptionName = "maxiter"; const std::string NativeEquationSolverSettings::maximalIterationsOptionShortName = "i"; const std::string NativeEquationSolverSettings::precisionOptionName = "precision"; const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute"; NativeEquationSolverSettings::NativeEquationSolverSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) { - std::vector<std::string> methods = { "jacobi" }; + std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor" }; this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine. Available are: { jacobi }.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("jacobi").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, omegaOptionName, false, "The omega used for SOR.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The value of the SOR parameter.").setDefaultValueDouble(0.9).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); } @@ -28,10 +31,14 @@ namespace storm { return this->getOption(techniqueOptionName).getHasOptionBeenSet(); } - NativeEquationSolverSettings::LinearEquationTechnique NativeEquationSolverSettings::getLinearEquationSystemTechnique() const { + NativeEquationSolverSettings::LinearEquationMethod NativeEquationSolverSettings::getLinearEquationSystemMethod() const { std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString(); if (linearEquationSystemTechniqueAsString == "jacobi") { - return NativeEquationSolverSettings::LinearEquationTechnique::Jacobi; + return NativeEquationSolverSettings::LinearEquationMethod::Jacobi; + } else if (linearEquationSystemTechniqueAsString == "gaussseidel") { + return NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel; + } else if (linearEquationSystemTechniqueAsString == "sor") { + return NativeEquationSolverSettings::LinearEquationMethod::SOR; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } @@ -51,6 +58,10 @@ namespace storm { double NativeEquationSolverSettings::getPrecision() const { return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble(); } + + double NativeEquationSolverSettings::getOmega() const { + return this->getOption(omegaOptionName).getArgumentByName("value").getValueAsDouble(); + } bool NativeEquationSolverSettings::isConvergenceCriterionSet() const { return this->getOption(absoluteOptionName).getHasOptionBeenSet(); diff --git a/src/settings/modules/NativeEquationSolverSettings.h b/src/settings/modules/NativeEquationSolverSettings.h index 9d6e43ecd..01c0fd7c9 100644 --- a/src/settings/modules/NativeEquationSolverSettings.h +++ b/src/settings/modules/NativeEquationSolverSettings.h @@ -12,8 +12,8 @@ namespace storm { */ class NativeEquationSolverSettings : public ModuleSettings { public: - // An enumeration of all available techniques for solving linear equations. - enum class LinearEquationTechnique { Jacobi }; + // An enumeration of all available methods for solving linear equations. + enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR }; // An enumeration of all available convergence criteria. enum class ConvergenceCriterion { Absolute, Relative }; @@ -33,11 +33,11 @@ namespace storm { bool isLinearEquationSystemTechniqueSet() const; /*! - * Retrieves the technique that is to be used for solving systems of linear equations. + * Retrieves the method that is to be used for solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - LinearEquationTechnique getLinearEquationSystemTechnique() const; + LinearEquationMethod getLinearEquationSystemMethod() const; /*! * Retrieves whether the maximal iteration count has been set. @@ -67,6 +67,13 @@ namespace storm { */ double getPrecision() const; + /*! + * Retrieves the value of omega to be used for the SOR method. + * + * @return The value of omega to be used. + */ + double getOmega() const; + /*! * Retrieves whether the convergence criterion has been set. * @@ -89,6 +96,7 @@ namespace storm { private: // Define the string names of the options as constants. static const std::string techniqueOptionName; + static const std::string omegaOptionName; static const std::string maximalIterationsOptionName; static const std::string maximalIterationsOptionShortName; static const std::string precisionOptionName; diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index 789ee6e55..bd43c3e88 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -31,31 +31,31 @@ namespace storm { restart = settings.getRestartIterationCount(); // Determine the method to be used. - storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique(); - if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab) { + storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod(); + if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab) { method = SolutionMethod::Bicgstab; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Qmr) { method = SolutionMethod::Qmr; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Gmres) { method = SolutionMethod::Gmres; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi) { method = SolutionMethod::Jacobi; } // Check which preconditioner to use. - storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique preconditionAsSetting = settings.getPreconditioningTechnique(); - if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu) { + storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod(); + if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Ilu) { preconditioner = Preconditioner::Ilu; - } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal) { + } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal) { preconditioner = Preconditioner::Diagonal; - } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::None) { + } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::None) { preconditioner = Preconditioner::None; } } template<typename ValueType> void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const { - LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'."); + LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations)."); if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) { LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored."); } diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index a2480ddc5..a389853c6 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -15,7 +15,7 @@ namespace storm { } template<typename ValueType> - NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) { + NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method) : A(A), method(method) { // Get the settings object to customize linear solving. storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); @@ -23,61 +23,94 @@ namespace storm { maximalNumberOfIterations = settings.getMaximalIterationCount(); precision = settings.getPrecision(); relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; - - // Determine the method to be used. - storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique(); - if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique::Jacobi) { - method = SolutionMethod::Jacobi; - } } template<typename ValueType> void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const { - // Get a Jacobi decomposition of the matrix A. - std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition(); - - // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. - bool multiplyResultProvided = true; - std::vector<ValueType>* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector<ValueType>(x.size()); - multiplyResultProvided = false; - } - std::vector<ValueType> const* copyX = nextX; - std::vector<ValueType>* currentX = &x; - - // Target vector for precision calculation. - std::vector<ValueType> tmpX(x.size()); - - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; - - while (!converged && iterationCount < maximalNumberOfIterations) { - // Compute D^-1 * (b - LU * x) and store result in nextX. - jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX); - storm::utility::vector::subtractVectors(b, tmpX, tmpX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX); + if (method == SolutionMethod::SOR || method == SolutionMethod::GaussSeidel) { + // Define the omega used for SOR. + ValueType omega = method == SolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one<ValueType>(); - // Swap the two pointers as a preparation for the next iteration. - std::swap(nextX, currentX); + // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. + bool tmpXProvided = true; + std::vector<ValueType>* tmpX = multiplyResult; + if (multiplyResult == nullptr) { + tmpX = new std::vector<ValueType>(x); + tmpXProvided = false; + } else { + *tmpX = x; + } - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative); + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } - - // If the last iteration did not write to the original x we have to swap the contents, because the - // output has to be written to the input parameter x. - if (currentX == copyX) { - std::swap(x, *currentX); - } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; + while (!converged && iterationCount < maximalNumberOfIterations) { + A.performSuccessiveOverRelaxationStep(omega, x, b); + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(precision), relative); + + // If we did not yet converge, we need to copy the contents of x to *tmpX. + if (!converged) { + *tmpX = x; + } + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (!tmpXProvided) { + delete tmpX; + } + } else { + // Get a Jacobi decomposition of the matrix A. + std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition(); + + // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. + bool multiplyResultProvided = true; + std::vector<ValueType>* nextX = multiplyResult; + if (nextX == nullptr) { + nextX = new std::vector<ValueType>(x.size()); + multiplyResultProvided = false; + } + std::vector<ValueType> const* copyX = nextX; + std::vector<ValueType>* currentX = &x; + + // Target vector for precision calculation. + std::vector<ValueType> tmpX(x.size()); + + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < maximalNumberOfIterations) { + // Compute D^-1 * (b - LU * x) and store result in nextX. + jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX); + storm::utility::vector::subtractVectors(b, tmpX, tmpX); + storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX); + + // Swap the two pointers as a preparation for the next iteration. + std::swap(nextX, currentX); + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative); + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // If the last iteration did not write to the original x we have to swap the contents, because the + // output has to be written to the input parameter x. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (!multiplyResultProvided) { + delete copyX; + } } } diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index f8f5a1378..aa20b1e71 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -14,15 +14,16 @@ namespace storm { public: // An enumeration specifying the available solution methods. enum class SolutionMethod { - Jacobi + Jacobi, GaussSeidel, SOR }; /*! * Constructs a linear equation solver with parameters being set according to the settings object. * * @param A The matrix defining the coefficients of the linear equation system. + * @param method The method to use for solving linear equations. */ - NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A); + NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method = SolutionMethod::Jacobi); /*! * Constructs a linear equation solver with the given parameters. diff --git a/src/solver/SymbolicGameSolver.cpp b/src/solver/SymbolicGameSolver.cpp new file mode 100644 index 000000000..cafa4eb04 --- /dev/null +++ b/src/solver/SymbolicGameSolver.cpp @@ -0,0 +1,69 @@ +#include "src/solver/SymbolicGameSolver.h" + +#include "src/storage/dd/CuddBdd.h" +#include "src/storage/dd/CuddAdd.h" + +namespace storm { + namespace solver { + + template<storm::dd::DdType Type> + SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables) { + // Get the settings object to customize solving. + storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); + storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings(); + + // Get appropriate settings. + maximalNumberOfIterations = settings.getMaximalIterationCount(); + precision = settings.getPrecision(); + relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; + } + + template<storm::dd::DdType Type> + SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) { + // Intentionally left empty. + } + + template<storm::dd::DdType Type> + storm::dd::Add<Type> SymbolicGameSolver<Type>::solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const { + // Set up the environment. + storm::dd::Add<Type> xCopy = x; + uint_fast64_t iterations = 0; + bool converged = false; + + while (!converged && iterations < maximalNumberOfIterations) { + // Compute tmp = A * x + b + storm::dd::Add<Type> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs); + storm::dd::Add<Type> tmp = this->gameMatrix.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables); + tmp += b; + + // Now abstract from player 2 and player 1 variables. + if (player2Min) { + tmp = tmp.minAbstract(this->player2Variables); + } else { + tmp = tmp.maxAbstract(this->player2Variables); + } + + if (player1Min) { + tmp = tmp.minAbstract(this->player1Variables); + } else { + tmp = tmp.maxAbstract(this->player1Variables); + } + + // Now check if the process already converged within our precision. + converged = xCopy.equalModuloPrecision(tmp, precision, relative); + + // If the method did not converge yet, we prepare the x vector for the next iteration. + if (!converged) { + xCopy = tmp; + } + + ++iterations; + } + + return xCopy; + } + + template class SymbolicGameSolver<storm::dd::DdType::CUDD>; + + } +} \ No newline at end of file diff --git a/src/solver/SymbolicGameSolver.h b/src/solver/SymbolicGameSolver.h new file mode 100644 index 000000000..ff5d64056 --- /dev/null +++ b/src/solver/SymbolicGameSolver.h @@ -0,0 +1,97 @@ +#ifndef STORM_SOLVER_SYMBOLICGAMESOLVER_H_ +#define STORM_SOLVER_SYMBOLICGAMESOLVER_H_ + +#include "src/storage/expressions/Variable.h" + +#include "src/storage/dd/Bdd.h" +#include "src/storage/dd/Add.h" + +namespace storm { + namespace solver { + + /*! + * An interface that represents an abstract symbolic game solver. + */ + template<storm::dd::DdType Type> + class SymbolicGameSolver { + public: + /*! + * Constructs a symbolic game solver with the given meta variable sets and pairs. + * + * @param gameMatrix The matrix defining the coefficients of the game. + * @param allRows A BDD characterizing all rows of the equation system. + * @param rowMetaVariables The meta variables used to encode the rows of the matrix. + * @param columnMetaVariables The meta variables used to encode the columns of the matrix. + * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta + * variables. + * @param player1Variables The meta variables used to encode the player 1 choices. + * @param player2Variables The meta variables used to encode the player 2 choices. + */ + SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables); + + /*! + * Constructs a symbolic game solver with the given meta variable sets and pairs. + * + * @param gameMatrix The matrix defining the coefficients of the game. + * @param allRows A BDD characterizing all rows of the equation system. + * @param rowMetaVariables The meta variables used to encode the rows of the matrix. + * @param columnMetaVariables The meta variables used to encode the columns of the matrix. + * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta + * variables. + * @param player1Variables The meta variables used to encode the player 1 choices. + * @param player2Variables The meta variables used to encode the player 2 choices. + * @param precision The precision to achieve. + * @param maximalNumberOfIterations The maximal number of iterations to perform when solving a linear + * equation system iteratively. + * @param relative Sets whether or not to use a relativ stopping criterion rather than an absolute one. + */ + SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables, double precision, uint_fast64_t maximalNumberOfIterations, bool relative); + + /*! + * Solves the equation system x = min/max(A*x + b) given by the parameters. Note that the matrix A has + * to be given upon construction time of the solver object. + * + * @param player1Min A flag indicating whether player 1 wants to minimize the result. + * @param player2Min A flag indicating whether player 1 wants to minimize the result. + * @param x The initial guess of the solution. + * @param b The vector to add after matrix-vector multiplication. + * @return The solution vector. + */ + virtual storm::dd::Add<Type> solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const; + + protected: + // The matrix defining the coefficients of the linear equation system. + storm::dd::Add<Type> const& gameMatrix; + + // A BDD characterizing all rows of the equation system. + storm::dd::Bdd<Type> const& allRows; + + // The row variables. + std::set<storm::expressions::Variable> rowMetaVariables; + + // The column variables. + std::set<storm::expressions::Variable> columnMetaVariables; + + // The pairs of meta variables used for renaming. + std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs; + + // The player 1 variables. + std::set<storm::expressions::Variable> player1Variables; + + // The player 2 variables. + std::set<storm::expressions::Variable> player2Variables; + + // The precision to achive. + double precision; + + // The maximal number of iterations to perform. + uint_fast64_t maximalNumberOfIterations; + + // A flag indicating whether a relative or an absolute stopping criterion is to be used. + bool relative; + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_SYMBOLICGAMESOLVER_H_ */ diff --git a/src/solver/SymbolicLinearEquationSolver.cpp b/src/solver/SymbolicLinearEquationSolver.cpp new file mode 100644 index 000000000..24a9d12b0 --- /dev/null +++ b/src/solver/SymbolicLinearEquationSolver.cpp @@ -0,0 +1,85 @@ +#include "src/solver/SymbolicLinearEquationSolver.h" + +#include "src/storage/dd/CuddDdManager.h" +#include "src/storage/dd/CuddAdd.h" + +namespace storm { + namespace solver { + + template<storm::dd::DdType DdType, typename ValueType> + SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) { + // Intentionally left empty. + } + + template<storm::dd::DdType DdType, typename ValueType> + SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) { + // Get the settings object to customize solving. + storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); + storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings(); + + // Get appropriate settings. + maximalNumberOfIterations = settings.getMaximalIterationCount(); + precision = settings.getPrecision(); + relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; + } + + template<storm::dd::DdType DdType, typename ValueType> + storm::dd::Add<DdType> SymbolicLinearEquationSolver<DdType, ValueType>::solveEquationSystem(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const { + // Start by computing the Jacobi decomposition of the matrix A. + storm::dd::Add<DdType> diagonal = x.getDdManager()->getAddOne(); + for (auto const& pair : rowColumnMetaVariablePairs) { + diagonal *= x.getDdManager()->getIdentity(pair.first).equals(x.getDdManager()->getIdentity(pair.second)); + diagonal *= x.getDdManager()->getRange(pair.first).toAdd() * x.getDdManager()->getRange(pair.second).toAdd(); + } + diagonal *= allRows.toAdd(); + + storm::dd::Add<DdType> lu = diagonal.ite(this->A.getDdManager()->getAddZero(), this->A); + storm::dd::Add<DdType> dinv = diagonal / (diagonal * this->A); + + // Set up additional environment variables. + storm::dd::Add<DdType> xCopy = x; + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < maximalNumberOfIterations) { + storm::dd::Add<DdType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs); + storm::dd::Add<DdType> tmp = lu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables); + tmp = b - tmp; + tmp = tmp.swapVariables(this->rowColumnMetaVariablePairs); + tmp = dinv.multiplyMatrix(tmp, this->columnMetaVariables); + + // Now check if the process already converged within our precision. + converged = xCopy.equalModuloPrecision(tmp, precision, relative); + + // If the method did not converge yet, we prepare the x vector for the next iteration. + if (!converged) { + xCopy = tmp; + } + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + return xCopy; + } + + template<storm::dd::DdType DdType, typename ValueType> + storm::dd::Add<DdType> SymbolicLinearEquationSolver<DdType, ValueType>::performMatrixVectorMultiplication(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b, uint_fast64_t n) const { + storm::dd::Add<DdType> xCopy = x; + + // Perform matrix-vector multiplication while the bound is met. + for (uint_fast64_t i = 0; i < n; ++i) { + xCopy = xCopy.swapVariables(this->rowColumnMetaVariablePairs); + xCopy = this->A.multiplyMatrix(xCopy, this->columnMetaVariables); + if (b != nullptr) { + xCopy += *b; + } + } + + return xCopy; + } + + template class SymbolicLinearEquationSolver<storm::dd::DdType::CUDD, double>; + + } +} diff --git a/src/solver/SymbolicLinearEquationSolver.h b/src/solver/SymbolicLinearEquationSolver.h new file mode 100644 index 000000000..6cd437673 --- /dev/null +++ b/src/solver/SymbolicLinearEquationSolver.h @@ -0,0 +1,104 @@ +#ifndef STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_ + +#include <memory> +#include <set> +#include <boost/variant.hpp> + +#include "src/storage/expressions/Variable.h" +#include "src/storage/dd/Bdd.h" +#include "src/storage/dd/Add.h" +#include "src/storage/dd/Odd.h" + +namespace storm { + namespace solver { + /*! + * An interface that represents an abstract symbolic linear equation solver. In addition to solving a system of + * linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided. + */ + template<storm::dd::DdType DdType, typename ValueType> + class SymbolicLinearEquationSolver { + public: + /*! + * Constructs a symbolic linear equation solver with the given meta variable sets and pairs. + * + * @param A The matrix defining the coefficients of the linear equation system. + * @param diagonal An ADD characterizing the elements on the diagonal of the matrix. + * @param allRows A BDD characterizing all rows of the equation system. + * @param rowMetaVariables The meta variables used to encode the rows of the matrix. + * @param columnMetaVariables The meta variables used to encode the columns of the matrix. + * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta + * variables. + */ + SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs); + + /*! + * Constructs a symbolic linear equation solver with the given meta variable sets and pairs. + * + * @param A The matrix defining the coefficients of the linear equation system. + * @param allRows A BDD characterizing all rows of the equation system. + * @param rowMetaVariables The meta variables used to encode the rows of the matrix. + * @param columnMetaVariables The meta variables used to encode the columns of the matrix. + * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta + * variables. + * @param precision The precision to achieve. + * @param maximalNumberOfIterations The maximal number of iterations to perform when solving a linear + * equation system iteratively. + * @param relative Sets whether or not to use a relativ stopping criterion rather than an absolute one. + */ + SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative); + + /*! + * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution. + * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has + * to be given upon construction time of the solver object. + * + * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A. + * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A. + * @return The solution of the equation system. + */ + virtual storm::dd::Add<DdType> solveEquationSystem(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const; + + /*! + * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After + * performing the necessary multiplications, the result is written to the input vector x. Note that the + * matrix A has to be given upon construction time of the solver object. + * + * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal + * to the number of rows of A. + * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal + * to the number of rows of A. + * @return The solution of the equation system. + */ + virtual storm::dd::Add<DdType> performMatrixVectorMultiplication(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b = nullptr, uint_fast64_t n = 1) const; + + protected: + // The matrix defining the coefficients of the linear equation system. + storm::dd::Add<DdType> const& A; + + // A BDD characterizing all rows of the equation system. + storm::dd::Bdd<DdType> const& allRows; + + // The row variables. + std::set<storm::expressions::Variable> rowMetaVariables; + + // The column variables. + std::set<storm::expressions::Variable> columnMetaVariables; + + // The pairs of meta variables used for renaming. + std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs; + + // The precision to achive. + double precision; + + // The maximal number of iterations to perform. + uint_fast64_t maximalNumberOfIterations; + + // A flag indicating whether a relative or an absolute stopping criterion is to be used. + bool relative; + }; + + } // namespace solver +} // namespace storm + +#endif /* STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_ */ diff --git a/src/solver/Z3SmtSolver.h b/src/solver/Z3SmtSolver.h index 1873cb4df..f00f67b89 100644 --- a/src/solver/Z3SmtSolver.h +++ b/src/solver/Z3SmtSolver.h @@ -75,7 +75,7 @@ namespace storm { * @param model The Z3 model to convert. * @return The valuation of variables corresponding to the given model. */ - storm::expressions::SimpleValuation convertZ3ModelToValuation(z3::model const& model); + storm::expressions::SimpleValuation convertZ3ModelToValuation(z3::model const& model); // The context used by the solver. std::unique_ptr<z3::context> context; @@ -87,10 +87,10 @@ namespace storm { std::unique_ptr<storm::adapters::Z3ExpressionAdapter> expressionAdapter; // A flag storing whether the last call to a check method provided aussumptions. - bool lastCheckAssumptions; + bool lastCheckAssumptions; // The last result that was returned by any of the check methods. - CheckResult lastResult; + CheckResult lastResult; #endif }; } diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 6684db2a8..ce5419099 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -50,7 +50,7 @@ namespace storm { void MatrixEntry<IndexType, ValueType>::setValue(ValueType const& value) { this->entry.second = value; } - + template<typename IndexType, typename ValueType> std::pair<IndexType, ValueType> const& MatrixEntry<IndexType, ValueType>::getColumnValuePair() const { return this->entry; @@ -82,13 +82,27 @@ namespace storm { rowIndications.push_back(0); } + template<typename ValueType> + SparseMatrixBuilder<ValueType>::SparseMatrixBuilder(SparseMatrix<ValueType>&& matrix) : initialRowCountSet(false), initialRowCount(0), initialColumnCountSet(false), initialColumnCount(0), initialEntryCountSet(false), initialEntryCount(0), forceInitialDimensions(false), hasCustomRowGrouping(matrix.nontrivialRowGrouping), initialRowGroupCountSet(false), initialRowGroupCount(0), rowGroupIndices(), columnsAndValues(std::move(matrix.columnsAndValues)), rowIndications(std::move(matrix.rowIndications)), currentEntryCount(matrix.entryCount), lastRow(matrix.rowCount - 1), lastColumn(columnsAndValues.back().getColumn()), highestColumn(matrix.getColumnCount() - 1), currentRowGroup() { + + // If the matrix has a custom row grouping, we move it and remove the last element to make it 'open' again. + if (hasCustomRowGrouping) { + rowGroupIndices = std::move(matrix.rowGroupIndices); + rowGroupIndices.pop_back(); + currentRowGroup = rowGroupIndices.size() - 1; + } + + // Likewise, we need to 'open' the row indications again. + rowIndications.pop_back(); + } + template<typename ValueType> void SparseMatrixBuilder<ValueType>::addNextValue(index_type row, index_type column, ValueType const& value) { // Check that we did not move backwards wrt. the row. if (row < lastRow) { throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in row " << row << ", but an element in row " << lastRow << " has already been added."; } - + // Check that we did not move backwards wrt. to column. if (row == lastRow && column < lastColumn) { throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in column " << column << " in row " << row << ", but an element in column " << lastColumn << " has already been added in that row."; @@ -103,14 +117,14 @@ namespace storm { lastRow = row; } - + lastColumn = column; // Finally, set the element and increase the current size. columnsAndValues.emplace_back(column, value); highestColumn = std::max(highestColumn, column); ++currentEntryCount; - + // In case we did not expect this value, we throw an exception. if (forceInitialDimensions) { STORM_LOG_THROW(!initialRowCountSet || lastRow < initialRowCount, storm::exceptions::OutOfRangeException, "Cannot insert value at illegal row " << lastRow << "."); @@ -145,14 +159,14 @@ namespace storm { // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for // the first and last row. rowIndications.push_back(currentEntryCount); - + uint_fast64_t columnCount = highestColumn + 1; if (initialColumnCountSet && forceInitialDimensions) { STORM_LOG_THROW(columnCount <= initialColumnCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialColumnCount << " columns, but got " << columnCount << "."); columnCount = std::max(columnCount, initialColumnCount); } columnCount = std::max(columnCount, overriddenColumnCount); - + uint_fast64_t entryCount = currentEntryCount; if (initialEntryCountSet && forceInitialDimensions) { STORM_LOG_THROW(entryCount == initialEntryCount, storm::exceptions::InvalidStateException, "Expected " << initialEntryCount << " entries, but got " << entryCount << "."); @@ -170,13 +184,23 @@ namespace storm { rowGroupCount = std::max(rowGroupCount, initialRowGroupCount); } rowGroupCount = std::max(rowGroupCount, overriddenRowGroupCount); - + for (index_type i = currentRowGroup; i <= rowGroupCount; ++i) { rowGroupIndices.push_back(rowCount); } } - - return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); + + return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), hasCustomRowGrouping); + } + + template<typename ValueType> + typename SparseMatrixBuilder<ValueType>::index_type SparseMatrixBuilder<ValueType>::getLastRow() const { + return lastRow; + } + + template<typename ValueType> + typename SparseMatrixBuilder<ValueType>::index_type SparseMatrixBuilder<ValueType>::getLastColumn() const { + return lastColumn; } template<typename ValueType> @@ -220,17 +244,24 @@ namespace storm { } template<typename ValueType> - SparseMatrix<ValueType>::SparseMatrix() : rowCount(0), columnCount(0), entryCount(0), nonzeroEntryCount(0), columnsAndValues(), rowIndications(), rowGroupIndices() { + SparseMatrix<ValueType>::SparseMatrix() : rowCount(0), columnCount(0), entryCount(0), nonzeroEntryCount(0), columnsAndValues(), rowIndications(), nontrivialRowGrouping(false), rowGroupIndices() { // Intentionally left empty. } template<typename ValueType> - SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType> const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), rowGroupIndices(other.rowGroupIndices) { + SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType> const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(other.rowGroupIndices) { // Intentionally left empty. } template<typename ValueType> - SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType>&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), rowGroupIndices(std::move(other.rowGroupIndices)) { + SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements) { + storm::storage::BitVector rowConstraint(other.getRowCount(), true); + storm::storage::BitVector columnConstraint(other.getColumnCount(), true); + *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements); + } + + template<typename ValueType> + SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType>&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(std::move(other.rowGroupIndices)) { // Now update the source matrix other.rowCount = 0; other.columnCount = 0; @@ -238,15 +269,15 @@ namespace storm { } template<typename ValueType> - SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) { - this->updateNonzeroEntryCount(); + SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) { + this->updateNonzeroEntryCount(); } template<typename ValueType> - SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) { - this->updateNonzeroEntryCount(); + SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(std::move(rowGroupIndices)) { + this->updateNonzeroEntryCount(); } - + template<typename ValueType> SparseMatrix<ValueType>& SparseMatrix<ValueType>::operator=(SparseMatrix<ValueType> const& other) { // Only perform assignment if source and target are not the same. @@ -259,6 +290,7 @@ namespace storm { columnsAndValues = other.columnsAndValues; rowIndications = other.rowIndications; rowGroupIndices = other.rowGroupIndices; + nontrivialRowGrouping = other.nontrivialRowGrouping; } return *this; @@ -276,6 +308,7 @@ namespace storm { columnsAndValues = std::move(other.columnsAndValues); rowIndications = std::move(other.rowIndications); rowGroupIndices = std::move(other.rowGroupIndices); + nontrivialRowGrouping = other.nontrivialRowGrouping; } return *this; @@ -288,7 +321,7 @@ namespace storm { } bool equalityResult = true; - + equalityResult &= this->getRowCount() == other.getRowCount(); if (!equalityResult) { return false; @@ -346,35 +379,35 @@ namespace storm { return entryCount; } - template<typename T> - uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const { - uint_fast64_t result = 0; - for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) { - result += (this->rowIndications[row + 1] - this->rowIndications[row]); - } - return result; - } - + template<typename T> + uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const { + uint_fast64_t result = 0; + for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) { + result += (this->rowIndications[row + 1] - this->rowIndications[row]); + } + return result; + } + template<typename ValueType> typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getNonzeroEntryCount() const { return nonzeroEntryCount; - } - - 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>()) { - ++this->nonzeroEntryCount; - } - } - } - - template<typename ValueType> - void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) { - this->nonzeroEntryCount += difference; - } + } + + 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>()) { + ++this->nonzeroEntryCount; + } + } + } + + template<typename ValueType> + void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) { + this->nonzeroEntryCount += difference; + } template<typename ValueType> typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getRowGroupCount() const { @@ -390,7 +423,7 @@ namespace storm { std::vector<typename SparseMatrix<ValueType>::index_type> const& SparseMatrix<ValueType>::getRowGroupIndices() const { return rowGroupIndices; } - + template<typename ValueType> void SparseMatrix<ValueType>::makeRowsAbsorbing(storm::storage::BitVector const& rows) { for (auto row : rows) { @@ -424,7 +457,7 @@ namespace storm { columnValuePtr->setValue(storm::utility::one<ValueType>()); ++columnValuePtr; for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) { - ++this->nonzeroEntryCount; + ++this->nonzeroEntryCount; columnValuePtr->setColumn(0); columnValuePtr->setValue(storm::utility::zero<ValueType>()); } @@ -477,12 +510,47 @@ namespace storm { return getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements); } } - + template<typename ValueType> SparseMatrix<ValueType> SparseMatrix<ValueType>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<index_type> const& rowGroupIndices, bool insertDiagonalEntries) const { - // First, we need to determine the number of entries and the number of rows of the submatrix. + uint_fast64_t submatrixColumnCount = columnConstraint.getNumberOfSetBits(); + + // 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>* 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; + } + } + + // Then, we need to determine the number of entries and the number of rows of the submatrix. index_type subEntries = 0; index_type subRows = 0; + index_type rowGroupCount = 0; for (auto index : rowGroupConstraint) { subRows += rowGroupIndices[index + 1] - rowGroupIndices[index]; for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { @@ -492,69 +560,55 @@ namespace storm { if (columnConstraint.get(it->getColumn())) { ++subEntries; - if (index == it->getColumn()) { + if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) { foundDiagonalElement = true; } } } // If requested, we need to reserve one entry more for inserting the diagonal zero entry. - if (insertDiagonalEntries && !foundDiagonalElement) { + if (insertDiagonalEntries && !foundDiagonalElement && rowGroupCount < submatrixColumnCount) { ++subEntries; } } + ++rowGroupCount; } // Create and initialize resulting matrix. - SparseMatrixBuilder<ValueType> matrixBuilder(subRows, columnConstraint.getNumberOfSetBits(), subEntries, true, true); - - // Create 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> bitsSetBeforeIndex; - bitsSetBeforeIndex.reserve(columnCount); - - // Compute the information to fill this vector. - index_type lastIndex = 0; - index_type currentNumberOfSetBits = 0; - - // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows are also - // taken. - storm::storage::BitVector columnBitCountConstraint = columnConstraint; - if (insertDiagonalEntries) { - columnBitCountConstraint |= rowGroupConstraint; - } - for (auto index : columnBitCountConstraint) { - while (lastIndex <= index) { - bitsSetBeforeIndex.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } + SparseMatrixBuilder<ValueType> matrixBuilder(subRows, submatrixColumnCount, subEntries, true, this->hasNontrivialRowGrouping()); // Copy over selected entries. + rowGroupCount = 0; index_type rowCount = 0; for (auto index : rowGroupConstraint) { - matrixBuilder.newRowGroup(rowCount); + if (this->hasNontrivialRowGrouping()) { + matrixBuilder.newRowGroup(rowCount); + } for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { bool insertedDiagonalElement = false; for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) { if (columnConstraint.get(it->getColumn())) { - if (index == it->getColumn()) { + if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) { insertedDiagonalElement = true; - } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) { - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>()); + } else if (insertDiagonalEntries && !insertedDiagonalElement && columnBitsSetBeforeIndex[it->getColumn()] > (*rowBitsSetBeforeIndex)[index]) { + matrixBuilder.addNextValue(rowCount, rowGroupCount, storm::utility::zero<ValueType>()); insertedDiagonalElement = true; } - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue()); + matrixBuilder.addNextValue(rowCount, columnBitsSetBeforeIndex[it->getColumn()], it->getValue()); } } - if (insertDiagonalEntries && !insertedDiagonalElement) { - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>()); + if (insertDiagonalEntries && !insertedDiagonalElement && rowGroupCount < submatrixColumnCount) { + matrixBuilder.addNextValue(rowGroupCount, rowGroupCount, storm::utility::zero<ValueType>()); } - ++rowCount; } + ++rowGroupCount; + } + + // If the constraints were not the same, we have allocated some additional memory we need to free now. + if (&rowGroupConstraint != &columnConstraint) { + delete rowBitsSetBeforeIndex; } return matrixBuilder.build(); @@ -615,12 +669,12 @@ namespace storm { SparseMatrix<ValueType> SparseMatrix<ValueType>::transpose(bool joinGroups, bool keepZeros) const { index_type rowCount = this->getColumnCount(); index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount(); - if (keepZeros) { - index_type entryCount = this->getEntryCount(); - } else { - this->updateNonzeroEntryCount(); - index_type entryCount = this->getNonzeroEntryCount(); - } + if (keepZeros) { + index_type entryCount = this->getEntryCount(); + } else { + this->updateNonzeroEntryCount(); + index_type entryCount = this->getNonzeroEntryCount(); + } std::vector<index_type> rowIndications(rowCount + 1); std::vector<MatrixEntry<index_type, ValueType>> columnsAndValues(entryCount); @@ -628,7 +682,7 @@ namespace storm { // First, we need to count how many entries each column has. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) { + if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) { ++rowIndications[transition.getColumn() + 1]; } } @@ -647,7 +701,7 @@ namespace storm { // Now we are ready to actually fill in the values of the transposed matrix. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) { + if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) { columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); nextIndices[transition.getColumn()]++; } @@ -659,11 +713,11 @@ namespace storm { rowGroupIndices[i] = i; } - storm::storage::SparseMatrix<ValueType> transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); - + storm::storage::SparseMatrix<ValueType> transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), false); + return transposedMatrix; } - + template<typename ValueType> void SparseMatrix<ValueType>::convertToEquationSystem() { invertDiagonal(); @@ -673,22 +727,22 @@ namespace storm { template<typename ValueType> void SparseMatrix<ValueType>::invertDiagonal() { // Now iterate over all row groups and set the diagonal elements to the inverted value. - // If there is a row without the diagonal element, an exception is thrown. - ValueType one = storm::utility::one<ValueType>(); - ValueType zero = storm::utility::zero<ValueType>(); + // If there is a row without the diagonal element, an exception is thrown. + ValueType one = storm::utility::one<ValueType>(); + ValueType zero = storm::utility::zero<ValueType>(); bool foundDiagonalElement = false; for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { - if (entry.getValue() == one) { - --this->nonzeroEntryCount; - entry.setValue(zero); - } else if (entry.getValue() == zero) { - ++this->nonzeroEntryCount; - entry.setValue(one); - } else { - entry.setValue(one - entry.getValue()); - } + if (entry.getValue() == one) { + --this->nonzeroEntryCount; + entry.setValue(zero); + } else if (entry.getValue() == zero) { + ++this->nonzeroEntryCount; + entry.setValue(one); + } else { + entry.setValue(one - entry.getValue()); + } foundDiagonalElement = true; } } @@ -718,7 +772,7 @@ namespace storm { for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { - --this->nonzeroEntryCount; + --this->nonzeroEntryCount; entry.setValue(storm::utility::zero<ValueType>()); } } @@ -799,15 +853,28 @@ namespace storm { typename std::vector<ValueType>::iterator resultIterator = result.begin(); typename std::vector<ValueType>::iterator resultIteratorEnd = result.end(); - for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { - *resultIterator = storm::utility::zero<ValueType>(); - - for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { - *resultIterator += it->getValue() * vector[it->getColumn()]; + // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable + // to store the intermediate result. + if (&vector == &result) { + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { + ValueType tmpValue = storm::utility::zero<ValueType>(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + tmpValue += it->getValue() * vector[it->getColumn()]; + } + *resultIterator = tmpValue; + } + } else { + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { + *resultIterator = storm::utility::zero<ValueType>(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + *resultIterator += it->getValue() * vector[it->getColumn()]; + } } } } - + #ifdef STORM_HAVE_INTELTBB template<typename ValueType> void SparseMatrix<ValueType>::multiplyWithVectorParallel(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const { @@ -832,7 +899,53 @@ namespace storm { }); } #endif - + + template<typename ValueType> + void SparseMatrix<ValueType>::performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const { + const_iterator it = this->begin(); + const_iterator ite; + std::vector<index_type>::const_iterator rowIterator = rowIndications.begin(); + typename std::vector<ValueType>::const_iterator bIt = b.begin(); + typename std::vector<ValueType>::const_iterator bIte = b.end(); + typename std::vector<ValueType>::iterator resultIterator = x.begin(); + typename std::vector<ValueType>::iterator resultIteratorEnd = x.end(); + + // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable + // to store the intermediate result. + index_type currentRow = 0; + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++bIt) { + ValueType tmpValue = storm::utility::zero<ValueType>(); + ValueType diagonalElement = storm::utility::zero<ValueType>(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + if (it->getColumn() != currentRow) { + tmpValue += it->getValue() * x[it->getColumn()]; + } else { + diagonalElement += it->getValue(); + } + } + + *resultIterator = ((storm::utility::one<ValueType>() - omega) * *resultIterator) + (omega / diagonalElement) * (*bIt - tmpValue); + ++currentRow; + } + } + + template<typename ValueType> + void SparseMatrix<ValueType>::multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const { + const_iterator it = this->begin(); + const_iterator ite; + std::vector<index_type>::const_iterator rowIterator = rowIndications.begin(); + std::vector<index_type>::const_iterator rowIteratorEnd = rowIndications.end(); + + uint_fast64_t currentRow = 0; + for (; rowIterator != rowIteratorEnd - 1; ++rowIterator) { + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + result[it->getColumn()] += it->getValue() * vector[currentRow]; + } + ++currentRow; + } + } + template<typename ValueType> std::size_t SparseMatrix<ValueType>::getSizeInBytes() const { uint_fast64_t size = sizeof(*this); @@ -906,6 +1019,11 @@ namespace storm { return this->columnsAndValues.begin() + this->rowIndications[rowCount]; } + template<typename ValueType> + bool SparseMatrix<ValueType>::hasNontrivialRowGrouping() const { + return nontrivialRowGrouping; + } + template<typename ValueType> ValueType SparseMatrix<ValueType>::getRowSum(index_type row) const { ValueType sum = storm::utility::zero<ValueType>(); @@ -921,7 +1039,7 @@ namespace storm { if (this->getRowCount() != matrix.getRowCount()) return false; if (this->getColumnCount() != matrix.getColumnCount()) return false; if (this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false; - + // Check the subset property for all rows individually. for (index_type row = 0; row < this->getRowCount(); ++row) { for (const_iterator it1 = this->begin(row), ite1 = this->end(row), it2 = matrix.begin(row), ite2 = matrix.end(row); it1 != ite1; ++it1) { diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 45e7b19ae..df682d756 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -149,6 +149,14 @@ namespace storm { */ SparseMatrixBuilder(index_type rows = 0, index_type columns = 0, index_type entries = 0, bool forceDimensions = true, bool hasCustomRowGrouping = false, index_type rowGroups = 0); + /*! + * Moves the contents of the given matrix into the matrix builder so that its contents can be modified again. + * This is, for example, useful if rows need to be added to the matrix. + * + * @param matrix The matrix that is to be made editable again. + */ + SparseMatrixBuilder(SparseMatrix<ValueType>&& matrix); + /*! * Sets the matrix entry at the given row and column to the given value. After all entries have been added, * a call to finalize(false) is mandatory. @@ -193,6 +201,20 @@ namespace storm { */ SparseMatrix<value_type> build(index_type overriddenRowCount = 0, index_type overriddenColumnCount = 0, index_type overriddenRowGroupCount = 0); + /*! + * Retrieves the most recently used row. + * + * @return The most recently used row. + */ + index_type getLastRow() const; + + /*! + * Retrieves the most recently used row. + * + * @return The most recently used row. + */ + index_type getLastColumn() const; + private: // A flag indicating whether a row count was set upon construction. bool initialRowCountSet; @@ -277,6 +299,7 @@ namespace storm { friend class storm::adapters::EigenAdapter; friend class storm::adapters::StormAdapter; friend class storm::solver::TopologicalValueIterationMinMaxLinearEquationSolver<ValueType>; + friend class SparseMatrixBuilder<ValueType>; typedef uint_fast64_t index_type; typedef ValueType value_type; @@ -388,6 +411,15 @@ namespace storm { */ SparseMatrix(SparseMatrix<value_type> const& other); + /*! + * Constructs a sparse matrix by performing a deep-copy of the given matrix. + * + * @param other The matrix from which to copy the content. + * @param insertDiagonalElements If set to true, the copy will have all diagonal elements. If they did not + * exist in the original matrix, they are inserted and set to value zero. + */ + SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements); + /*! * Constructs a sparse matrix by moving the contents of the given matrix to the newly created one. * @@ -402,8 +434,9 @@ namespace storm { * @param rowIndications The row indications vector of the matrix to be constructed. * @param columnsAndValues The vector containing the columns and values of the entries in the matrix. * @param rowGroupIndices The vector representing the row groups in the matrix. + * @param hasNontrivialRowGrouping If set to true, this indicates that the row grouping is non-trivial. */ - SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices); + SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool hasNontrivialRowGrouping); /*! * Constructs a sparse matrix by moving the given contents. @@ -412,8 +445,9 @@ namespace storm { * @param rowIndications The row indications vector of the matrix to be constructed. * @param columnsAndValues The vector containing the columns and values of the entries in the matrix. * @param rowGroupIndices The vector representing the row groups in the matrix. + * @param hasNontrivialRowGrouping If set to true, this indicates that the row grouping is non-trivial. */ - SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, value_type>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices); + SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, value_type>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices, bool hasNontrivialRowGrouping); /*! * Assigns the contents of the given matrix to the current one by deep-copying its contents. @@ -646,6 +680,26 @@ namespace storm { * @return The product of the matrix and the given vector as the content of the given result vector. */ void multiplyWithVector(std::vector<value_type> const& vector, std::vector<value_type>& result) const; + + /*! + * Multiplies the vector to the matrix from the left and writes the result to the given result vector. + * + * @param vector The vector with which the matrix is to be multiplied. This vector is interpreted as being + * a row vector. + * @param result The vector that is supposed to hold the result of the multiplication after the operation. + * @return The product of the matrix and the given vector as the content of the given result vector. The + * result is to be interpreted as a row vector. + */ + void multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const; + + /*! + * Performs one step of the successive over-relaxation technique. + * + * @param omega The Omega parameter for SOR. + * @param x The current solution vector. The result will be written to the very same vector. + * @param b The 'right-hand side' of the problem. + */ + void performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const; /*! * Multiplies the matrix with the given vector in a sequential way and writes the result to the given result @@ -798,6 +852,13 @@ namespace storm { * @return An iterator that points past the end of the last row of the matrix. */ iterator end(); + + /*! + * Retrieves whether the matrix has a (possibly) non-trivial row grouping. + * + * @return True iff the matrix has a (possibly) non-trivial row grouping. + */ + bool hasNontrivialRowGrouping() const; /*! * Returns a copy of the matrix with the chosen internal data type @@ -813,7 +874,7 @@ namespace storm { new_columnsAndValues.at(i) = MatrixEntry<SparseMatrix::index_type, NewValueType>(columnsAndValues.at(i).getColumn(), static_cast<NewValueType>(columnsAndValues.at(i).getValue())); } - return SparseMatrix<NewValueType>(columnCount, std::move(new_rowIndications), std::move(new_columnsAndValues), std::move(new_rowGroupIndices)); + return SparseMatrix<NewValueType>(columnCount, std::move(new_rowIndications), std::move(new_columnsAndValues), std::move(new_rowGroupIndices), nontrivialRowGrouping); } private: @@ -852,6 +913,10 @@ namespace storm { // entry is not included anymore. std::vector<index_type> rowIndications; + // A flag that indicates whether the matrix has a non-trivial row-grouping, i.e. (possibly) more than one + // row per row group. + bool nontrivialRowGrouping; + // A vector indicating the row groups of the matrix. std::vector<index_type> rowGroupIndices; }; diff --git a/src/storage/dd/CuddAdd.cpp b/src/storage/dd/CuddAdd.cpp index 6b7f2067f..5eaded696 100644 --- a/src/storage/dd/CuddAdd.cpp +++ b/src/storage/dd/CuddAdd.cpp @@ -557,7 +557,7 @@ namespace storm { rowIndications[0] = 0; // Construct matrix and return result. - return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(trivialRowGroupIndices)); + return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(trivialRowGroupIndices), false); } storm::storage::SparseMatrix<double> Add<DdType::CUDD>::toMatrix(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const { @@ -680,7 +680,7 @@ namespace storm { } rowIndications[0] = 0; - return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); + return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), true); } std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> Add<DdType::CUDD>::toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const { @@ -804,7 +804,7 @@ namespace storm { } rowIndications[0] = 0; - return std::make_pair(storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)), std::move(explicitVector)); + return std::make_pair(storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), true), std::move(explicitVector)); } template<typename ValueType> diff --git a/src/storage/dd/CuddAdd.h b/src/storage/dd/CuddAdd.h index 01a57b629..03c8b3a2c 100644 --- a/src/storage/dd/CuddAdd.h +++ b/src/storage/dd/CuddAdd.h @@ -721,7 +721,7 @@ namespace storm { void toMatrixRec(DdNode const* dd, std::vector<uint_fast64_t>& rowIndications, std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>>& columnsAndValues, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, Odd<DdType::CUDD> const& columnOdd, uint_fast64_t currentRowLevel, uint_fast64_t currentColumnLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, uint_fast64_t currentColumnOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices, std::vector<uint_fast64_t> const& ddColumnVariableIndices, bool generateValues = true) const; /*! - * Helper function to convert the DD into a (sparse) vector. + * Helper function to convert the DD into an (explicit) vector. * * @param dd The DD to convert. * @param result The vector that will hold the values upon successful completion. @@ -772,7 +772,7 @@ namespace storm { void addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const; /*! - * Performs a recursive step to perform the given function between the given DD-based vector to the given + * Performs a recursive step to perform the given function between the given DD-based vector and the given * explicit vector. * * @param dd The DD to add to the explicit vector. diff --git a/src/storage/dd/CuddBdd.cpp b/src/storage/dd/CuddBdd.cpp index ca93c1742..30ed8a191 100644 --- a/src/storage/dd/CuddBdd.cpp +++ b/src/storage/dd/CuddBdd.cpp @@ -337,6 +337,41 @@ namespace storm { } } + storm::storage::BitVector Bdd<DdType::CUDD>::toVector(storm::dd::Odd<DdType::CUDD> const& rowOdd) const { + std::vector<uint_fast64_t> ddVariableIndices = this->getSortedVariableIndices(); + storm::storage::BitVector result(rowOdd.getTotalOffset()); + this->toVectorRec(this->getCuddDdNode(), this->getDdManager()->getCuddManager(), result, rowOdd, Cudd_IsComplement(this->getCuddDdNode()), 0, ddVariableIndices.size(), 0, ddVariableIndices); + return result; + } + + void Bdd<DdType::CUDD>::toVectorRec(DdNode const* dd, Cudd const& manager, storm::storage::BitVector& result, Odd<DdType::CUDD> const& rowOdd, bool complement, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const { + // If there are no more values to select, we can directly return. + if (dd == Cudd_ReadLogicZero(manager.getManager()) && !complement) { + return; + } else if (dd == Cudd_ReadOne(manager.getManager()) && complement) { + return; + } + + // If we are at the maximal level, the value to be set is stored as a constant in the DD. + if (currentRowLevel == maxLevel) { + result.set(currentRowOffset, true); + } else if (ddRowVariableIndices[currentRowLevel] < dd->index) { + toVectorRec(dd, manager, result, rowOdd.getElseSuccessor(), complement, currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices); + toVectorRec(dd, manager, result, rowOdd.getThenSuccessor(), complement, currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices); + } else { + // Otherwise, we compute the ODDs for both the then- and else successors. + DdNode* elseDdNode = Cudd_E(dd); + DdNode* thenDdNode = Cudd_T(dd); + + // Determine whether we have to evaluate the successors as if they were complemented. + bool elseComplemented = Cudd_IsComplement(elseDdNode) ^ complement; + bool thenComplemented = Cudd_IsComplement(thenDdNode) ^ complement; + + toVectorRec(Cudd_Regular(elseDdNode), manager, result, rowOdd.getElseSuccessor(), elseComplemented, currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices); + toVectorRec(Cudd_Regular(thenDdNode), manager, result, rowOdd.getThenSuccessor(), thenComplemented, currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices); + } + } + std::ostream& operator<<(std::ostream& out, const Bdd<DdType::CUDD>& bdd) { bdd.exportToDot(); return out; diff --git a/src/storage/dd/CuddBdd.h b/src/storage/dd/CuddBdd.h index 326072eb9..ad2d9aa6a 100644 --- a/src/storage/dd/CuddBdd.h +++ b/src/storage/dd/CuddBdd.h @@ -263,6 +263,15 @@ namespace storm { */ Add<DdType::CUDD> toAdd() const; + /*! + * Converts the BDD to a bit vector. The given offset-labeled DD is used to determine the correct row of + * each entry. + * + * @param rowOdd The ODD used for determining the correct row. + * @return The bit vector that is represented by this BDD. + */ + storm::storage::BitVector toVector(storm::dd::Odd<DdType::CUDD> const& rowOdd) const; + private: /*! * Retrieves the CUDD BDD object associated with this DD. @@ -315,6 +324,21 @@ namespace storm { template<typename ValueType> static DdNode* fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::function<bool (ValueType const&)> const& filter); + /*! + * Helper function to convert the DD into a bit vector. + * + * @param dd The DD to convert. + * @param manager The Cudd manager responsible for the DDs. + * @param result The vector that will hold the values upon successful completion. + * @param rowOdd The ODD used for the row translation. + * @param complement A flag indicating whether the result is to be interpreted as a complement. + * @param currentRowLevel The currently considered row level in the DD. + * @param maxLevel The number of levels that need to be considered. + * @param currentRowOffset The current row offset. + * @param ddRowVariableIndices The (sorted) indices of all DD row variables that need to be considered. + */ + void toVectorRec(DdNode const* dd, Cudd const& manager, storm::storage::BitVector& result, Odd<DdType::CUDD> const& rowOdd, bool complement, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const; + // The BDD created by CUDD. BDD cuddBdd; }; diff --git a/src/storage/expressions/ExprtkExpressionEvaluator.cpp b/src/storage/expressions/ExprtkExpressionEvaluator.cpp index b621b64a6..4e79a8dc0 100644 --- a/src/storage/expressions/ExprtkExpressionEvaluator.cpp +++ b/src/storage/expressions/ExprtkExpressionEvaluator.cpp @@ -49,7 +49,7 @@ namespace storm { CompiledExpressionType& compiledExpression = result.first->second; compiledExpression.register_symbol_table(symbolTable); bool parsingOk = parser.compile(ToExprtkStringVisitor().toString(expression), compiledExpression); - STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk."); + STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk: " << *expression); return compiledExpression; } diff --git a/src/utility/cli.cpp b/src/utility/cli.cpp index 50c15e627..7e7a99080 100644 --- a/src/utility/cli.cpp +++ b/src/utility/cli.cpp @@ -34,11 +34,11 @@ namespace storm { } void printHeader(const int argc, const char* argv[]) { - std::cout << "StoRM" << std::endl; + std::cout << "StoRM" << std::endl; std::cout << "--------" << std::endl << std::endl; - // std::cout << storm::utility::StormVersion::longVersionString() << std::endl; + std::cout << storm::utility::StormVersion::longVersionString() << std::endl; #ifdef STORM_HAVE_INTELTBB std::cout << "Linked with Intel Threading Building Blocks v" << TBB_VERSION_MAJOR << "." << TBB_VERSION_MINOR << " (Interface version " << TBB_INTERFACE_VERSION << ")." << std::endl; #endif @@ -113,6 +113,7 @@ namespace storm { std::cout << "Command line arguments: " << commandStream.str() << std::endl; std::cout << "Current working directory: " << getCurrentWorkingDirectory() << std::endl << std::endl; } + } void printUsage() { diff --git a/src/utility/cli.h b/src/utility/cli.h index 5b924ea78..82044200b 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -251,6 +251,22 @@ namespace storm { } #ifdef STORM_HAVE_CARL + void exportParametricResultToFile(storm::RationalFunction const& result, storm::models::sparse::Dtmc<storm::RationalFunction>::ConstraintCollector const& constraintCollector, std::string const& path) { + std::ofstream filestream; + filestream.open(path); + // TODO: add checks. + filestream << "!Parameters: "; + std::set<storm::Variable> vars = result.gatherVariables(); + std::copy(vars.begin(), vars.end(), std::ostream_iterator<storm::Variable>(filestream, ", ")); + filestream << std::endl; + filestream << "!Result: " << result << std::endl; + filestream << "!Well-formed Constraints: " << std::endl; + std::copy(constraintCollector.getWellformedConstraints().begin(), constraintCollector.getWellformedConstraints().end(), std::ostream_iterator<carl::Constraint<storm::RationalFunction>>(filestream, "\n")); + filestream << "!Graph-preserving Constraints: " << std::endl; + std::copy(constraintCollector.getGraphPreservingConstraints().begin(), constraintCollector.getGraphPreservingConstraints().end(), std::ostream_iterator<carl::Constraint<storm::RationalFunction>>(filestream, "\n")); + filestream.close(); + } + template<> inline void verifySparseModel(boost::optional<storm::prism::Program> const& program, std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>> model, std::shared_ptr<storm::logic::Formula> formula) { storm::settings::modules::GeneralSettings const& settings = storm::settings::generalSettings(); @@ -276,6 +292,11 @@ namespace storm { } else { std::cout << " skipped, because the modelling formalism is currently unsupported." << std::endl; } + + storm::settings::modules::ParametricSettings const& parametricSettings = storm::settings::parametricSettings(); + if (parametricSettings.exportResultToFile()) { + exportParametricResultToFile(result->asExplicitQuantitativeCheckResult<storm::RationalFunction>()[*dtmc->getInitialStates().begin()], storm::models::sparse::Dtmc<storm::RationalFunction>::ConstraintCollector(*dtmc), parametricSettings.exportResultPath()); + } } #endif @@ -462,4 +483,4 @@ namespace storm { } } -#endif \ No newline at end of file +#endif diff --git a/src/utility/graph.h b/src/utility/graph.h index 4cb349a4e..45d90dffb 100644 --- a/src/utility/graph.h +++ b/src/utility/graph.h @@ -245,7 +245,7 @@ namespace storm { * with probability 0 and the second stores all indices of states with probability 1. */ template <typename T> - static std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01(storm::storage::SparseMatrix<T> backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) { + static std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01(storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) { std::pair<storm::storage::BitVector, storm::storage::BitVector> result; result.first = performProbGreater0(backwardTransitions, phiStates, psiStates); result.second = performProb1(backwardTransitions, phiStates, psiStates, result.first); diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp index 2d2c0cdcb..785f1d00c 100644 --- a/src/utility/solver.cpp +++ b/src/utility/solver.cpp @@ -2,6 +2,8 @@ #include "src/settings/SettingsManager.h" +#include "src/solver/SymbolicGameSolver.h" + #include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/GmmxxLinearEquationSolver.h" @@ -15,6 +17,16 @@ namespace storm { namespace utility { namespace solver { + template<storm::dd::DdType Type, typename ValueType> + std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> SymbolicLinearEquationSolverFactory<Type, ValueType>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const { + return std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicLinearEquationSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs)); + } + + template<storm::dd::DdType Type> + std::unique_ptr<storm::solver::SymbolicGameSolver<Type>> SymbolicGameSolverFactory<Type>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const { + return std::unique_ptr<storm::solver::SymbolicGameSolver<Type>>(new storm::solver::SymbolicGameSolver<Type>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables)); + } + template<typename ValueType> std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const { storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver(); @@ -29,9 +41,28 @@ namespace storm { return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix)); } + template<typename ValueType> + NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory() { + switch (storm::settings::nativeEquationSolverSettings().getLinearEquationSystemMethod()) { + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi: + this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::Jacobi; + break; + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel: + this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::GaussSeidel; + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR: + this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::SOR; + } + omega = storm::settings::nativeEquationSolverSettings().getOmega(); + } + + template<typename ValueType> + NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega) : method(method), omega(omega) { + // Intentionally left empty. + } + template<typename ValueType> std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const { - return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix)); + return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix, method)); } template<typename ValueType> @@ -79,6 +110,8 @@ namespace storm { return factory->create(name); } + template class SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>; + template class SymbolicGameSolverFactory<storm::dd::DdType::CUDD>; template class LinearEquationSolverFactory<double>; template class GmmxxLinearEquationSolverFactory<double>; template class NativeLinearEquationSolverFactory<double>; diff --git a/src/utility/solver.h b/src/utility/solver.h index a72bbdbb1..cd56695f0 100644 --- a/src/utility/solver.h +++ b/src/utility/solver.h @@ -1,15 +1,34 @@ #ifndef STORM_UTILITY_SOLVER_H_ #define STORM_UTILITY_SOLVER_H_ +#include "src/solver/SymbolicGameSolver.h" + +#include "src/solver/SymbolicLinearEquationSolver.h" #include "src/solver/LinearEquationSolver.h" +#include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/MinMaxLinearEquationSolver.h" #include "src/solver/LpSolver.h" +#include "src/storage/dd/DdType.h" +#include "src/settings/modules/NativeEquationSolverSettings.h" + #include "src/exceptions/InvalidSettingsException.h" namespace storm { namespace utility { namespace solver { + template<storm::dd::DdType Type, typename ValueType> + class SymbolicLinearEquationSolverFactory { + public: + virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const; + }; + + template<storm::dd::DdType Type> + class SymbolicGameSolverFactory { + public: + virtual std::unique_ptr<storm::solver::SymbolicGameSolver<Type>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const; + }; + template<typename ValueType> class LinearEquationSolverFactory { public: @@ -25,7 +44,14 @@ namespace storm { template<typename ValueType> class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> { public: + NativeLinearEquationSolverFactory(); + NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega); + virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override; + + private: + typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method; + ValueType omega; }; template<typename ValueType> diff --git a/src/utility/vector.h b/src/utility/vector.h index b5567ae5b..634d2bc5a 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -376,7 +376,9 @@ namespace storm { if (val2 == 0) { return (std::abs(val1) <= precision); } - if (std::abs((val1 - val2)/val2) > precision) return false; + if (std::abs((val1 - val2)/val2) > precision) { + return false; + } } else { if (std::abs(val1 - val2) > precision) return false; } diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp index 4c79aa929..80300f1f4 100644 --- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp @@ -86,6 +86,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>(); EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"minimum\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult8 = checkResult->asExplicitQuantitativeCheckResult<double>(); + EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8[initialState], storm::settings::generalSettings().getPrecision()); } TEST(GmmxxCtmcCslModelCheckerTest, Embedded) { @@ -149,6 +156,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>(); EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"fail_sensors\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>(); + EXPECT_NEAR(0.93458866427696596, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision()); } TEST(GmmxxCtmcCslModelCheckerTest, Polling) { @@ -176,6 +190,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Polling) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult1 = checkResult->asExplicitQuantitativeCheckResult<double>(); EXPECT_NEAR(1, quantitativeCheckResult1[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=?[\"target\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult2 = checkResult->asExplicitQuantitativeCheckResult<double>(); + EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2[initialState], storm::settings::generalSettings().getPrecision()); } TEST(GmmxxCtmcCslModelCheckerTest, Fms) { @@ -252,4 +273,11 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>(); EXPECT_NEAR(262.85103661561755, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"first_queue_full\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>(); + EXPECT_NEAR(0.9100373532, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision()); } diff --git a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp index 08fb5a29c..9bbedfa7b 100644 --- a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp @@ -127,3 +127,148 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) { EXPECT_NEAR(1.044879046, quantitativeResult3[0], storm::settings::gmmxxEquationSolverSettings().getPrecision()); } + +TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) { + storm::storage::SparseMatrixBuilder<double> matrixBuilder; + std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2); + matrixBuilder.addNextValue(0, 1, 1.); + matrixBuilder.addNextValue(1, 0, 1.); + storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); + auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + { + matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4); + matrixBuilder.addNextValue(0, 0, .5); + matrixBuilder.addNextValue(0, 1, .5); + matrixBuilder.addNextValue(1, 0, .5); + matrixBuilder.addNextValue(1, 1, .5); + storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); + auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + + { + matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 2, 1); + matrixBuilder.addNextValue(2, 0, 1); + storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 2); + + dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); + auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>(); + + EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + } +} + +TEST(GmmxxDtmcPrctlModelCheckerTest, LRA) { + storm::storage::SparseMatrixBuilder<double> matrixBuilder; + std::shared_ptr<storm::models::sparse::Dtmc<double>> mdp; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder<double>(15, 15, 20, true); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 4, 0.7); + matrixBuilder.addNextValue(1, 6, 0.3); + matrixBuilder.addNextValue(2, 0, 1); + + matrixBuilder.addNextValue(3, 5, 0.8); + matrixBuilder.addNextValue(3, 9, 0.2); + matrixBuilder.addNextValue(4, 3, 1); + matrixBuilder.addNextValue(5, 3, 1); + + matrixBuilder.addNextValue(6, 7, 1); + matrixBuilder.addNextValue(7, 8, 1); + matrixBuilder.addNextValue(8, 6, 1); + + matrixBuilder.addNextValue(9, 10, 1); + matrixBuilder.addNextValue(10, 9, 1); + matrixBuilder.addNextValue(11, 9, 1); + + matrixBuilder.addNextValue(12, 5, 0.4); + matrixBuilder.addNextValue(12, 8, 0.3); + matrixBuilder.addNextValue(12, 11, 0.3); + + matrixBuilder.addNextValue(13, 7, 0.7); + matrixBuilder.addNextValue(13, 12, 0.3); + + matrixBuilder.addNextValue(14, 12, 1); + + storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(15); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + ap.addLabelToState("a", 4); + ap.addLabelToState("a", 5); + ap.addLabelToState("a", 7); + ap.addLabelToState("a", 11); + ap.addLabelToState("a", 13); + ap.addLabelToState("a", 14); + + mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); + auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>(); + + EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3/3., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision()); + } +} diff --git a/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp index 968a9476b..a720507a4 100644 --- a/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp @@ -102,6 +102,14 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Cluster) { storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"minimum\"]"); + checkResult = modelchecker.check(*formula); + + checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates())); + storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult8 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8.getMin(), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8.getMax(), storm::settings::generalSettings().getPrecision()); } TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded) { @@ -173,6 +181,14 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded) { storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"fail_sensors\"]"); + checkResult = modelchecker.check(*formula); + + checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates())); + storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + EXPECT_NEAR(0.934586179, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.934586179, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision()); } TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling) { @@ -201,6 +217,15 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling) { storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); EXPECT_NEAR(1, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(1, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"target\"]"); + checkResult = modelchecker.check(*formula); + + checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates())); + storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.20079750055570736, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision()); + } TEST(GmmxxHybridCtmcCslModelCheckerTest, Fms) { @@ -288,4 +313,12 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Tandem) { storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"first_queue_full\"]"); + checkResult = modelchecker.check(*formula); + + checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates())); + storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + EXPECT_NEAR(0.9100373532, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.9100373532, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision()); } diff --git a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp index 7f53f45b4..f58b4259e 100644 --- a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "storm-config.h" +#include "src/settings/SettingMemento.h" #include "src/logic/Formulas.h" #include "src/utility/solver.h" #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" @@ -128,7 +129,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, SynchronousLeader) { EXPECT_NEAR(1.0448979589010925, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision()); } -TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) { +TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) { storm::storage::SparseMatrixBuilder<double> matrixBuilder; std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc; @@ -144,7 +145,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); @@ -169,7 +170,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); @@ -194,7 +195,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); @@ -255,7 +256,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA) { mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>())); + storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a"); auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula); diff --git a/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp new file mode 100644 index 000000000..02d6ad5ba --- /dev/null +++ b/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp @@ -0,0 +1,172 @@ +#include "gtest/gtest.h" +#include "storm-config.h" + +#include "src/logic/Formulas.h" +#include "src/utility/solver.h" +#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h" +#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h" +#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h" +#include "src/parser/PrismParser.h" +#include "src/builder/DdPrismModelBuilder.h" +#include "src/models/symbolic/Dtmc.h" +#include "src/settings/SettingsManager.h" + +TEST(SymbolicDtmcPrctlModelCheckerTest, Die) { + storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/die.pm"); + + // Build the die model with its reward model. +#ifdef WINDOWS + storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options; +#else + typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options; +#endif + options.buildRewards = true; + options.rewardModelName = "coin_flips"; + std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options); + EXPECT_EQ(13, model->getNumberOfStates()); + EXPECT_EQ(20, model->getNumberOfTransitions()); + + ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc); + + std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>(); + + storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("one"); + auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two"); + eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three"); + eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + auto done = std::make_shared<storm::logic::AtomicLabelFormula>("done"); + auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(done); + + result = checker.check(*reachabilityRewardFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(3.6666622161865234, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(3.6666622161865234, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); +} + +TEST(SymbolicDtmcPrctlModelCheckerTest, Crowds) { + storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/crowds-5-5.pm"); + + std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program); + EXPECT_EQ(8607, model->getNumberOfStates()); + EXPECT_EQ(15113, model->getNumberOfTransitions()); + + ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc); + + std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>(); + + storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observe0Greater1"); + auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(0.33288236360191303, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.33288236360191303, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeIGreater1"); + eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(0.15222081144084315, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.15222081144084315, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeOnlyTrueSender"); + eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(0.3215392962289586, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3215392962289586, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); +} + +TEST(SymbolicDtmcPrctlModelCheckerTest, SynchronousLeader) { + storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader-3-5.pm"); + + // Build the die model with its reward model. +#ifdef WINDOWS + storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options; +#else + typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options; +#endif + options.buildRewards = true; + options.rewardModelName = "num_rounds"; + std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options); + EXPECT_EQ(273, model->getNumberOfStates()); + EXPECT_EQ(397, model->getNumberOfTransitions()); + + ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc); + + std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>(); + + storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>())); + + auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected"); + auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula); + + std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(1.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected"); + auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true); + auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 20); + + result = checker.check(*boundedUntilFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected"); + auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula); + + result = checker.check(*reachabilityRewardFormula); + result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates())); + storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>(); + + EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision()); +} + diff --git a/test/functional/solver/FullySymbolicGameSolverTest.cpp b/test/functional/solver/FullySymbolicGameSolverTest.cpp new file mode 100644 index 000000000..e41606d7f --- /dev/null +++ b/test/functional/solver/FullySymbolicGameSolverTest.cpp @@ -0,0 +1,64 @@ +#include "gtest/gtest.h" +#include "storm-config.h" + +#include "src/storage/dd/CuddDdManager.h" + +#include "src/utility/solver.h" + +TEST(FullySymbolicGameSolverTest, Solve) { + // Create some variables. + std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::CUDD>> manager(new storm::dd::DdManager<storm::dd::DdType::CUDD>()); + std::pair<storm::expressions::Variable, storm::expressions::Variable> state = manager->addMetaVariable("x", 1, 4); + std::pair<storm::expressions::Variable, storm::expressions::Variable> pl1 = manager->addMetaVariable("a", 0, 1); + std::pair<storm::expressions::Variable, storm::expressions::Variable> pl2 = manager->addMetaVariable("b", 0, 1); + + storm::dd::Bdd<storm::dd::DdType::CUDD> allRows = manager->getBddZero(); + std::set<storm::expressions::Variable> rowMetaVariables({state.first}); + std::set<storm::expressions::Variable> columnMetaVariables({state.second}); + std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> rowColumnMetaVariablePairs = {state}; + std::set<storm::expressions::Variable> player1Variables({pl1.first}); + std::set<storm::expressions::Variable> player2Variables({pl2.first}); + + // Construct simple game. + storm::dd::Add<storm::dd::DdType::CUDD> matrix = manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 2).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.6); + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 1).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.4); + + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 2).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(0.2); + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 3).toAdd() * manager->getEncoding(pl1.first, 0).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(0.8); + + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 3).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.5); + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 4).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 0).toAdd() * manager->getConstant(0.5); + + matrix += manager->getEncoding(state.first, 1).toAdd() * manager->getEncoding(state.second, 1).toAdd() * manager->getEncoding(pl1.first, 1).toAdd() * manager->getEncoding(pl2.first, 1).toAdd() * manager->getConstant(1); + + std::unique_ptr<storm::utility::solver::SymbolicGameSolverFactory<storm::dd::DdType::CUDD>> solverFactory(new storm::utility::solver::SymbolicGameSolverFactory<storm::dd::DdType::CUDD>()); + std::unique_ptr<storm::solver::SymbolicGameSolver<storm::dd::DdType::CUDD>> solver = solverFactory->create(matrix, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables,player2Variables); + + // Create solution and target state vector. + storm::dd::Add<storm::dd::DdType::CUDD> x = manager->getAddZero(); + storm::dd::Add<storm::dd::DdType::CUDD> b = manager->getEncoding(state.first, 2).toAdd() + manager->getEncoding(state.first, 4).toAdd(); + + // Now solve the game with different strategies for the players. + storm::dd::Add<storm::dd::DdType::CUDD> result = solver->solveGame(true, true, x, b); + result *= manager->getEncoding(state.first, 1).toAdd(); + result = result.sumAbstract({state.first}); + EXPECT_NEAR(0, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + x = manager->getAddZero(); + result = solver->solveGame(true, false, x, b); + result *= manager->getEncoding(state.first, 1).toAdd(); + result = result.sumAbstract({state.first}); + EXPECT_NEAR(0.5, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + x = manager->getAddZero(); + result = solver->solveGame(false, true, x, b); + result *= manager->getEncoding(state.first, 1).toAdd(); + result = result.sumAbstract({state.first}); + EXPECT_NEAR(0.2, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision()); + + x = manager->getAddZero(); + result = solver->solveGame(false, false, x, b); + result *= manager->getEncoding(state.first, 1).toAdd(); + result = result.sumAbstract({state.first}); + EXPECT_NEAR(0.99999892625817599, result.getValue(), storm::settings::nativeEquationSolverSettings().getPrecision()); +} \ No newline at end of file diff --git a/test/functional/storage/SparseMatrixTest.cpp b/test/functional/storage/SparseMatrixTest.cpp index 20716916a..d16b5558d 100644 --- a/test/functional/storage/SparseMatrixTest.cpp +++ b/test/functional/storage/SparseMatrixTest.cpp @@ -154,8 +154,8 @@ TEST(SparseMatrix, CreationWithMovingContents) { columnsAndValues.emplace_back(1, 0.7); columnsAndValues.emplace_back(3, 0.2); - ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3})); - storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}); + ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}, false)); + storm::storage::SparseMatrix<double> matrix(4, {0, 2, 5, 5}, columnsAndValues, {0, 1, 2, 3}, false); ASSERT_EQ(3, matrix.getRowCount()); ASSERT_EQ(4, matrix.getColumnCount()); ASSERT_EQ(5, matrix.getEntryCount()); @@ -328,7 +328,7 @@ TEST(SparseMatrix, Submatrix) { ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 3, 0.3)); storm::storage::SparseMatrix<double> matrix; ASSERT_NO_THROW(matrix = matrixBuilder.build()); - + std::vector<uint_fast64_t> rowGroupIndices = {0, 1, 2, 4, 5}; storm::storage::BitVector rowGroupConstraint(4);