From 7e672cddd9b1a5e613f16ff4a165dc7a8026bcda Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Tue, 17 Mar 2015 17:50:21 +0100 Subject: [PATCH 01/20] Started implementation of LRA for MDPs - adapted storm::utility::graph::getReachableStates to work for non-deterministic matrices Former-commit-id: cd7e4697575e60eea6ee5e5dc9361cf599207c46 --- .../prctl/SparseMdpPrctlModelChecker.cpp | 151 ++++++++++++++++++ .../prctl/SparseMdpPrctlModelChecker.h | 4 + src/utility/graph.h | 2 +- 3 files changed, 156 insertions(+), 1 deletion(-) diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index ced5463a6..18c4f268e 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -12,6 +12,9 @@ #include "src/exceptions/InvalidPropertyException.h" +#include "src/storage/MaximalEndComponentDecomposition.h" +#include "src/storage/MaximalEndComponent.h" + namespace storm { namespace modelchecker { template @@ -313,6 +316,154 @@ namespace storm { return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); } + + template + std::vector SparseMdpPrctlModelChecker::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const { + // If there are no goal states, we avoid the computation and directly return zero. + if (psiStates.empty()) { + return std::vector(model.getNumberOfStates(), storm::utility::zero()); + } + + // Likewise, if all bits are set, we can avoid the computation and set. + if ((~psiStates).empty()) { + return std::vector(model.getNumberOfStates(), storm::utility::one()); + } + + // Start by decomposing the MDP into its MECs. + storm::storage::MaximalEndComponentDecomposition mecDecomposition(model); + + // Get some data members for convenience. + typename storm::storage::SparseMatrix const& transitionMatrix = model.getTransitionMatrix(); + std::vector const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); + + // Now start with compute the long-run average for all end components in isolation. + std::vector lraValuesForEndComponents; + + // While doing so, we already gather some information for the following steps. + std::vector stateToMecIndexMap(model.getNumberOfStates()); + storm::storage::BitVector statesInMecs(model.getNumberOfStates()); + + for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { + storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; + + // Gather information for later use. + for (auto const& stateChoicesPair : mec) { + uint_fast64_t state = stateChoicesPair.first; + + statesInMecs.set(state); + stateToMecIndexMap[state] = currentMecIndex; + } + + // Compute the LRA value for the current MEC. + lraValuesForEndComponents.push_back(this->computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec)); + } + + // For fast transition rewriting, we build some auxiliary data structures. + storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; + + // Prepare result vector. + std::vector result(model.getNumberOfStates()); + + //Set the values for all states in MECs. + for (auto state : statesInMecs) { + result[state] = lraValuesForEndComponents[stateToMecIndexMap[state]]; + } + + //for all states not in any mec set the result to the minimal/maximal value of the reachable MECs + //there might be a more efficient way to do this... + for (auto state : statesNotContainedInAnyMec){ + + //calculate what result values the reachable states in MECs have + storm::storage::BitVector currentState(model.getNumberOfStates); + currentState.set(state); + storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( + transitionMatrix, currentState, storm::storage::BitVector(model.getNumberOfStates, true), statesInMecs + ); + + storm::storage::BitVector reachableMecStates = statesInMecs & reachableStates; + std::vector reachableResults(reachableMecStates.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(reachableResults, reachableMecStates, result); + + //now select the minimal/maximal element + if (minimize){ + result[state] = *std::min_element(reachableResults.begin(), reachableResults.end()); + } else { + result[state] = *std::max_element(reachableResults.begin(), reachableResults.end()); + } + } + + return result; + } + + template + std::unique_ptr SparseMdpPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); + + std::unique_ptr subResultPointer = this->check(stateFormula); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); + } + + template + ValueType SparseMarkovAutomatonCslModelChecker::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { + std::shared_ptr solver = storm::utility::solver::getLpSolver("LRA for MEC"); + solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize); + + //// First, we need to create the variables for the problem. + //std::map stateToVariableMap; + //for (auto const& stateChoicesPair : mec) { + // std::string variableName = "x" + std::to_string(stateChoicesPair.first); + // stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName); + //} + //storm::expressions::Variable k = solver->addUnboundedContinuousVariable("k", 1); + //solver->update(); + + //// Now we encode the problem as constraints. + //for (auto const& stateChoicesPair : mec) { + // uint_fast64_t state = stateChoicesPair.first; + + // // Now, based on the type of the state, create a suitable constraint. + // if (markovianStates.get(state)) { + // storm::expressions::Expression constraint = stateToVariableMap.at(state); + + // for (auto element : transitionMatrix.getRow(nondeterministicChoiceIndices[state])) { + // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); + // } + + // constraint = constraint + solver->getConstant(storm::utility::one() / exitRates[state]) * k; + // storm::expressions::Expression rightHandSide = psiStates.get(state) ? solver->getConstant(storm::utility::one() / exitRates[state]) : solver->getConstant(0); + // if (minimize) { + // constraint = constraint <= rightHandSide; + // } else { + // constraint = constraint >= rightHandSide; + // } + // solver->addConstraint("state" + std::to_string(state), constraint); + // } else { + // // For probabilistic states, we want to add the constraint x_s <= sum P(s, a, s') * x_s' where a is the current action + // // and the sum ranges over all states s'. + // for (auto choice : stateChoicesPair.second) { + // storm::expressions::Expression constraint = stateToVariableMap.at(state); + + // for (auto element : transitionMatrix.getRow(choice)) { + // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); + // } + + // storm::expressions::Expression rightHandSide = solver->getConstant(storm::utility::zero()); + // if (minimize) { + // constraint = constraint <= rightHandSide; + // } else { + // constraint = constraint >= rightHandSide; + // } + // solver->addConstraint("state" + std::to_string(state), constraint); + // } + // } + //} + + //solver->optimize(); + //return solver->getContinuousValue(k); + } + template storm::models::Mdp const& SparseMdpPrctlModelChecker::getModel() const { return this->template getModelAs>(); diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 630ce0bb2..3dad8e619 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -39,6 +39,7 @@ namespace storm { virtual std::unique_ptr computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; protected: storm::models::Mdp const& getModel() const override; @@ -52,7 +53,10 @@ namespace storm { std::vector computeInstantaneousRewardsHelper(bool minimize, uint_fast64_t stepCount) const; std::vector computeCumulativeRewardsHelper(bool minimize, uint_fast64_t stepBound) const; std::vector computeReachabilityRewardsHelper(bool minimize, storm::storage::BitVector const& targetStates, bool qualitative) const; + std::vector computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const; + static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec); + // A solver that is used for solving systems of linear equations that are the result of nondeterministic choices. std::shared_ptr> nondeterministicLinearEquationSolver; }; diff --git a/src/utility/graph.h b/src/utility/graph.h index 95b9e4ddc..dc8a2c98c 100644 --- a/src/utility/graph.h +++ b/src/utility/graph.h @@ -51,7 +51,7 @@ namespace storm { currentState = stack.back(); stack.pop_back(); - for (auto const& successor : transitionMatrix.getRow(currentState)) { + for (auto const& successor : transitionMatrix.getRowGroup(currentState)) { // Only explore the state if the transition was actually there and the successor has not yet // been visited. if (successor.getValue() != storm::utility::zero() && !reachableStates.get(successor.getColumn())) { From 25739720e0219d0b57e073a5f335c690d3373acb Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 19 Mar 2015 15:33:41 +0100 Subject: [PATCH 02/20] Finished implementation of LRA for MPDs. No tests yet. Former-commit-id: 795c0e98425e63af91fe03ac5e5f3341ad70c16a --- .../prctl/SparseMdpPrctlModelChecker.cpp | 114 ++++++++---------- .../prctl/SparseMdpPrctlModelChecker.h | 2 + 2 files changed, 53 insertions(+), 63 deletions(-) diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 18c4f268e..cccc2bbe2 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -1,6 +1,7 @@ #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" #include +#include #include "src/utility/constants.h" #include "src/utility/macros.h" @@ -10,10 +11,12 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/solver/LpSolver.h" + #include "src/exceptions/InvalidPropertyException.h" +#include "src/storage/expressions/Expressions.h" #include "src/storage/MaximalEndComponentDecomposition.h" -#include "src/storage/MaximalEndComponent.h" namespace storm { namespace modelchecker { @@ -320,28 +323,28 @@ namespace storm { template std::vector SparseMdpPrctlModelChecker::computeLongRunAverageHelper(bool minimize, 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(model.getNumberOfStates(), storm::utility::zero()); + return std::vector(numOfStates, storm::utility::zero()); } // Likewise, if all bits are set, we can avoid the computation and set. if ((~psiStates).empty()) { - return std::vector(model.getNumberOfStates(), storm::utility::one()); + return std::vector(numOfStates, storm::utility::one()); } // Start by decomposing the MDP into its MECs. - storm::storage::MaximalEndComponentDecomposition mecDecomposition(model); + storm::storage::MaximalEndComponentDecomposition mecDecomposition(this->getModelAs>()); // Get some data members for convenience. - typename storm::storage::SparseMatrix const& transitionMatrix = model.getTransitionMatrix(); - std::vector const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); + typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); // Now start with compute the long-run average for all end components in isolation. std::vector lraValuesForEndComponents; // While doing so, we already gather some information for the following steps. - std::vector stateToMecIndexMap(model.getNumberOfStates()); - storm::storage::BitVector statesInMecs(model.getNumberOfStates()); + std::vector stateToMecIndexMap(numOfStates); + storm::storage::BitVector statesInMecs(numOfStates); for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; @@ -362,7 +365,7 @@ namespace storm { storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; // Prepare result vector. - std::vector result(model.getNumberOfStates()); + std::vector result(numOfStates); //Set the values for all states in MECs. for (auto state : statesInMecs) { @@ -374,10 +377,10 @@ namespace storm { for (auto state : statesNotContainedInAnyMec){ //calculate what result values the reachable states in MECs have - storm::storage::BitVector currentState(model.getNumberOfStates); + storm::storage::BitVector currentState(numOfStates); currentState.set(state); storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( - transitionMatrix, currentState, storm::storage::BitVector(model.getNumberOfStates, true), statesInMecs + transitionMatrix, currentState, storm::storage::BitVector(numOfStates, true), statesInMecs ); storm::storage::BitVector reachableMecStates = statesInMecs & reachableStates; @@ -406,62 +409,47 @@ namespace storm { } template - ValueType SparseMarkovAutomatonCslModelChecker::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { + ValueType SparseMdpPrctlModelChecker::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { std::shared_ptr solver = storm::utility::solver::getLpSolver("LRA for MEC"); solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize); //// First, we need to create the variables for the problem. - //std::map stateToVariableMap; - //for (auto const& stateChoicesPair : mec) { - // std::string variableName = "x" + std::to_string(stateChoicesPair.first); - // stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName); - //} - //storm::expressions::Variable k = solver->addUnboundedContinuousVariable("k", 1); - //solver->update(); - - //// Now we encode the problem as constraints. - //for (auto const& stateChoicesPair : mec) { - // uint_fast64_t state = stateChoicesPair.first; - - // // Now, based on the type of the state, create a suitable constraint. - // if (markovianStates.get(state)) { - // storm::expressions::Expression constraint = stateToVariableMap.at(state); - - // for (auto element : transitionMatrix.getRow(nondeterministicChoiceIndices[state])) { - // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); - // } - - // constraint = constraint + solver->getConstant(storm::utility::one() / exitRates[state]) * k; - // storm::expressions::Expression rightHandSide = psiStates.get(state) ? solver->getConstant(storm::utility::one() / exitRates[state]) : solver->getConstant(0); - // if (minimize) { - // constraint = constraint <= rightHandSide; - // } else { - // constraint = constraint >= rightHandSide; - // } - // solver->addConstraint("state" + std::to_string(state), constraint); - // } else { - // // For probabilistic states, we want to add the constraint x_s <= sum P(s, a, s') * x_s' where a is the current action - // // and the sum ranges over all states s'. - // for (auto choice : stateChoicesPair.second) { - // storm::expressions::Expression constraint = stateToVariableMap.at(state); - - // for (auto element : transitionMatrix.getRow(choice)) { - // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); - // } - - // storm::expressions::Expression rightHandSide = solver->getConstant(storm::utility::zero()); - // if (minimize) { - // constraint = constraint <= rightHandSide; - // } else { - // constraint = constraint >= rightHandSide; - // } - // solver->addConstraint("state" + std::to_string(state), constraint); - // } - // } - //} - - //solver->optimize(); - //return solver->getContinuousValue(k); + std::map stateToVariableMap; + for (auto const& stateChoicesPair : mec) { + std::string variableName = "h" + std::to_string(stateChoicesPair.first); + stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName); + } + storm::expressions::Variable lambda = solver->addUnboundedContinuousVariable("L", 1); + solver->update(); + + // Now we encode the problem as constraints. + for (auto const& stateChoicesPair : mec) { + uint_fast64_t state = stateChoicesPair.first; + + // Now, based on the type of the state, create a suitable constraint. + for (auto choice : stateChoicesPair.second) { + storm::expressions::Expression constraint = solver->getConstant(1); + ValueType w = 0; + + for (auto element : transitionMatrix.getRow(choice)) { + constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); + if (psiStates.get(element.getColumn())) { + w += element.getValue(); + } + } + constraint = constraint - solver->getConstant(w) * lambda; + + if (minimize) { + constraint = stateToVariableMap.at(state) <= constraint; + } else { + constraint = stateToVariableMap.at(state) >= constraint; + } + solver->addConstraint("state" + std::to_string(state) + "," + std::to_string(choice), constraint); + } + } + + solver->optimize(); + return solver->getContinuousValue(lambda); } template diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 3dad8e619..9718457cd 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -5,6 +5,8 @@ #include "src/models/Mdp.h" #include "src/utility/solver.h" #include "src/solver/NondeterministicLinearEquationSolver.h" +#include "src/storage/MaximalEndComponent.h" + namespace storm { namespace counterexamples { From b096180de808916f8dc2a7c845a01b8b015875da Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 20 Mar 2015 17:15:49 +0100 Subject: [PATCH 03/20] LRA on DTMCs implemented Former-commit-id: 633d81323d15dc040410cd73b682d4072bbf65fe --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 129 ++++++++++++++++++ .../prctl/SparseDtmcPrctlModelChecker.h | 7 +- 2 files changed, 135 insertions(+), 1 deletion(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 9d12d6166..c7aa308e0 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -1,6 +1,7 @@ #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include +#include #include "src/utility/macros.h" #include "src/utility/vector.h" @@ -9,6 +10,8 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/storage/StronglyConnectedComponentDecomposition.h" + #include "src/exceptions/InvalidPropertyException.h" namespace storm { @@ -297,6 +300,132 @@ namespace storm { return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeReachabilityRewardsHelper(subResult.getTruthValuesVector(), qualitative))); } + template + std::vector SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(bool minimize, 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(numOfStates, storm::utility::zero()); + } + + // Likewise, if all bits are set, we can avoid the computation and set. + if ((~psiStates).empty()) { + return std::vector(numOfStates, storm::utility::one()); + } + + // Start by decomposing the DTMC into its BSCCs. + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(this->getModelAs>(), false, true); + + // Get some data members for convenience. + typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); + + // Now start with compute the long-run average for all BSCCs in isolation. + std::vector lraValuesForBsccs; + + // While doing so, we already gather some information for the following steps. + std::vector stateToBsccIndexMap(numOfStates); + storm::storage::BitVector statesInBsccs(numOfStates); + + 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; + } + + // Compute the LRA value for the current BSCC + lraValuesForBsccs.push_back(this->computeLraForBSCC(transitionMatrix, psiStates, bscc)); + } + + // For fast transition rewriting, we build some auxiliary data structures. + storm::storage::BitVector statesNotContainedInAnyBscc = ~statesInBsccs; + + // Prepare result vector. + std::vector result(numOfStates); + + //Set the values for all states in BSCCs. + for (auto state : statesInBsccs) { + result[state] = lraValuesForBsccs[stateToBsccIndexMap[state]]; + } + + //for all states not in any bscc set the result to the minimal/maximal value of the reachable BSCCs + //there might be a more efficient way to do this... + for (auto state : statesNotContainedInAnyBscc){ + + //calculate what result values the reachable states in BSCCs have + storm::storage::BitVector currentState(numOfStates); + currentState.set(state); + storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( + transitionMatrix, currentState, storm::storage::BitVector(numOfStates, true), statesInBsccs + ); + + storm::storage::BitVector reachableBsccStates = statesInBsccs & reachableStates; + std::vector reachableResults(reachableBsccStates.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(reachableResults, reachableBsccStates, result); + + //now select the minimal/maximal element + if (minimize){ + result[state] = *std::min_element(reachableResults.begin(), reachableResults.end()); + } else { + result[state] = *std::max_element(reachableResults.begin(), reachableResults.end()); + } + } + + return result; + } + + template + std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); + + std::unique_ptr subResultPointer = this->check(stateFormula); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); + } + + + template + ValueType SparseDtmcPrctlModelChecker::computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc) { + //if size is one we can immediately derive the result + if (bscc.size() == 1){ + if (psiStates.get(*(bscc.begin()))) { + return storm::utility::one(); + } else{ + return storm::utility::zero(); + } + } + std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); + + storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount()); + subsystem.set(bscc.begin(), bscc.end()); + + storm::storage::SparseMatrix subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem); + std::vector steadyStateDist(subsystemMatrix.getRowCount(), storm::utility::zero()); + steadyStateDist[0] = storm::utility::one(); + + solver->solveEquationSystem(subsystemMatrix, steadyStateDist, steadyStateDist); + + ValueType length = storm::utility::zero(); + for (auto value : steadyStateDist) { + length += value; + } + storm::utility::vector::scaleVectorInPlace(steadyStateDist, storm::utility::one() / length); + + std::vector steadyStateForPsiStates(transitionMatrix.getRowCount(), storm::utility::zero()); + storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist); + + ValueType result = storm::utility::zero(); + + for (auto value : steadyStateForPsiStates) { + result += value; + } + + return result; + } + template storm::models::Dtmc const& SparseDtmcPrctlModelChecker::getModel() const { return this->template getModelAs>(); diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index f7557d63b..d1e880e1a 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -5,6 +5,7 @@ #include "src/models/Dtmc.h" #include "src/utility/solver.h" #include "src/solver/LinearEquationSolver.h" +#include "src/storage/StronglyConnectedComponent.h" namespace storm { namespace modelchecker { @@ -23,6 +24,7 @@ namespace storm { virtual std::unique_ptr computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; protected: storm::models::Dtmc const& getModel() const override; @@ -35,7 +37,10 @@ namespace storm { std::vector computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const; std::vector computeCumulativeRewardsHelper(uint_fast64_t stepBound) const; std::vector computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const; - + std::vector computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const; + + static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc); + // An object that is used for solving linear equations and performing matrix-vector multiplication. std::unique_ptr> linearEquationSolver; }; From 04c1d513137f4428f836b4f206130cff8a8e3b91 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Tue, 7 Apr 2015 13:59:39 +0200 Subject: [PATCH 04/20] intermediate commit, copied transpose and get submatrix code over and started adapting it. (changing workplace) Former-commit-id: af4a34dd3bb7ecdabdaac5cdec4b6353838c3e7a --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 123 +++++++++++++++++- 1 file changed, 122 insertions(+), 1 deletion(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index c7aa308e0..23d39c35d 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -402,6 +402,127 @@ namespace storm { 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. As a memory optimisation, we calculate + // ((P')^t - I) in one step from transitionMatrix instead of calling the appropriate methods + + // First, we need to determine the number of entries and the number of rows of the submatrix. + typename transitionMatrix::index_type subEntries = 0; + typename transitionMatrix::index_type subRows = 0; + for (auto index : rowGroupConstraint) { + subRows += rowGroupIndices[index + 1] - rowGroupIndices[index]; + for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { + bool foundDiagonalElement = false; + + for (const_iterator it = transitionMatrix->begin(i), ite = transitionMatrix->end(i); it != ite; ++it) { + if (columnConstraint.get(it->getColumn())) { + ++subEntries; + + if (index == it->getColumn()) { + foundDiagonalElement = true; + } + } + } + + // If requested, we need to reserve one entry more for inserting the diagonal zero entry. + if (!foundDiagonalElement) { + ++subEntries; + } + } + } + + // 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 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; + } + + // Copy over selected entries. + index_type rowCount = 0; + for (auto index : rowGroupConstraint) { + 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()) { + insertedDiagonalElement = true; + } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) { + matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero()); + insertedDiagonalElement = true; + } + matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue()); + } + } + if (insertDiagonalEntries && !insertedDiagonalElement) { + matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero()); + } + + ++rowCount; + } + } + matrixBuilder.build(); + + //transpose code + std::vector rowIndications(rowCount + 1); + std::vector> columnsAndValues(entryCount); + + // First, we need to count how many entries each column has. + for (index_type group = 0; group < columnCount; ++group) { + for (auto const& transition : transitionMatrix->getRow(group)) { + if (!comparator.isZero(transition.getValue())) { + ++rowIndications[transition.getColumn() + 1]; + } + } + } + + // Now compute the accumulated offsets. + for (index_type i = 1; i < rowCount + 1; ++i) { + rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; + } + + // Create an array that stores the index for the next value to be added for + // each row in the transposed matrix. Initially this corresponds to the previously + // computed accumulated offsets. + std::vector nextIndices = rowIndications; + + // 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 : transitionMatrix->getRow(group)) { + if (!comparator.isZero(transition.getValue())) { + columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); + nextIndices[transition.getColumn()]++; + } + } + } + + std::vector rowGroupIndices(rowCount + 1); + for (index_type i = 0; i <= rowCount; ++i) { + rowGroupIndices[i] = i; + } + + storm::storage::SparseMatrix transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); + + + /* storm::storage::SparseMatrix subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem); std::vector steadyStateDist(subsystemMatrix.getRowCount(), storm::utility::zero()); steadyStateDist[0] = storm::utility::one(); @@ -422,7 +543,7 @@ namespace storm { for (auto value : steadyStateForPsiStates) { result += value; } - + */ return result; } From a448cd8973b79b2010388eae81780e25b9513f4f Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 9 Apr 2015 21:57:00 +0200 Subject: [PATCH 05/20] Calculating steady state using standard equation system for eigenvectors, removed all-in-one matrix transformation (nicer looking code) Former-commit-id: 250261568636eea192f87f8a63e73739acae48eb --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 160 ++++-------------- 1 file changed, 34 insertions(+), 126 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 23d39c35d..4e35fef9c 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -403,147 +403,55 @@ namespace storm { 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. As a memory optimisation, we calculate - // ((P')^t - I) in one step from transitionMatrix instead of calling the appropriate methods + // ^t means transose, and I is the identity matrix. - // First, we need to determine the number of entries and the number of rows of the submatrix. - typename transitionMatrix::index_type subEntries = 0; - typename transitionMatrix::index_type subRows = 0; - for (auto index : rowGroupConstraint) { - subRows += rowGroupIndices[index + 1] - rowGroupIndices[index]; - for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { - bool foundDiagonalElement = false; - - for (const_iterator it = transitionMatrix->begin(i), ite = transitionMatrix->end(i); it != ite; ++it) { - if (columnConstraint.get(it->getColumn())) { - ++subEntries; - - if (index == it->getColumn()) { - foundDiagonalElement = true; - } - } - } - - // If requested, we need to reserve one entry more for inserting the diagonal zero entry. - if (!foundDiagonalElement) { - ++subEntries; + storm::storage::SparseMatrix 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 equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount()); + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); + 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()); } } - } - - // 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 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; + // 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."); } - - // Copy over selected entries. - index_type rowCount = 0; - for (auto index : rowGroupConstraint) { - 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()) { - insertedDiagonalElement = true; - } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) { - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero()); - insertedDiagonalElement = true; - } - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue()); - } - } - if (insertDiagonalEntries && !insertedDiagonalElement) { - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero()); - } - - ++rowCount; - } + for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) { + equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one); } - matrixBuilder.build(); + equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero); + subsystemMatrix = equationSystemBuilder.build(); - //transpose code - std::vector rowIndications(rowCount + 1); - std::vector> columnsAndValues(entryCount); + // 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 steadyStateDist(subsystemMatrix.getRowCount(), zero); + std::vector b(subsystemMatrix.getRowCount(), zero); + b[subsystemMatrix.getRowCount() - 1] = one; - // First, we need to count how many entries each column has. - for (index_type group = 0; group < columnCount; ++group) { - for (auto const& transition : transitionMatrix->getRow(group)) { - if (!comparator.isZero(transition.getValue())) { - ++rowIndications[transition.getColumn() + 1]; - } - } - } + solver->solveEquationSystem(subsystemMatrix, steadyStateDist, b); - // Now compute the accumulated offsets. - for (index_type i = 1; i < rowCount + 1; ++i) { - rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; - } - - // Create an array that stores the index for the next value to be added for - // each row in the transposed matrix. Initially this corresponds to the previously - // computed accumulated offsets. - std::vector nextIndices = rowIndications; - - // 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 : transitionMatrix->getRow(group)) { - if (!comparator.isZero(transition.getValue())) { - columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); - nextIndices[transition.getColumn()]++; - } - } - } - - std::vector rowGroupIndices(rowCount + 1); - for (index_type i = 0; i <= rowCount; ++i) { - rowGroupIndices[i] = i; - } - - storm::storage::SparseMatrix transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); - - - /* - storm::storage::SparseMatrix subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem); - std::vector steadyStateDist(subsystemMatrix.getRowCount(), storm::utility::zero()); - steadyStateDist[0] = storm::utility::one(); - - solver->solveEquationSystem(subsystemMatrix, steadyStateDist, steadyStateDist); - - ValueType length = storm::utility::zero(); - for (auto value : steadyStateDist) { - length += value; - } - storm::utility::vector::scaleVectorInPlace(steadyStateDist, storm::utility::one() / length); - - std::vector steadyStateForPsiStates(transitionMatrix.getRowCount(), storm::utility::zero()); + //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 steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero); storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist); - ValueType result = storm::utility::zero(); + ValueType result = zero; for (auto value : steadyStateForPsiStates) { result += value; } - */ + return result; } From 53f2fdf51e1387e1bc5e9d0d0c32a8ab3a6bbf65 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 10 Apr 2015 16:08:59 +0200 Subject: [PATCH 06/20] Changed implementation of LRA to be weighted with the probability to reach BSCCs instead of choosing min/max Former-commit-id: 347fda8e2226b08a30265b4ad39715454efc4889 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 66 +++++++------------ .../prctl/SparseDtmcPrctlModelChecker.h | 2 +- 2 files changed, 24 insertions(+), 44 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 4e35fef9c..1f8642b54 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -301,7 +301,7 @@ namespace storm { } template - std::vector SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const { + std::vector SparseDtmcPrctlModelChecker::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()) { @@ -319,71 +319,51 @@ namespace storm { // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); - // Now start with compute the long-run average for all BSCCs in isolation. - std::vector lraValuesForBsccs; - - // While doing so, we already gather some information for the following steps. - std::vector stateToBsccIndexMap(numOfStates); - storm::storage::BitVector statesInBsccs(numOfStates); + // First we check which states are in BSCCs. We use this later to speed up reachability analysis + storm::storage::BitVector statesNotInBsccs(numOfStates, true); 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; + statesNotInBsccs.set(state, false); } - - // Compute the LRA value for the current BSCC - lraValuesForBsccs.push_back(this->computeLraForBSCC(transitionMatrix, psiStates, bscc)); } - // For fast transition rewriting, we build some auxiliary data structures. - storm::storage::BitVector statesNotContainedInAnyBscc = ~statesInBsccs; - // Prepare result vector. std::vector result(numOfStates); - //Set the values for all states in BSCCs. - for (auto state : statesInBsccs) { - result[state] = lraValuesForBsccs[stateToBsccIndexMap[state]]; - } + //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]; - //for all states not in any bscc set the result to the minimal/maximal value of the reachable BSCCs - //there might be a more efficient way to do this... - for (auto state : statesNotContainedInAnyBscc){ - - //calculate what result values the reachable states in BSCCs have - storm::storage::BitVector currentState(numOfStates); - currentState.set(state); - storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( - transitionMatrix, currentState, storm::storage::BitVector(numOfStates, true), statesInBsccs - ); - - storm::storage::BitVector reachableBsccStates = statesInBsccs & reachableStates; - std::vector reachableResults(reachableBsccStates.getNumberOfSetBits()); - storm::utility::vector::selectVectorValues(reachableResults, reachableBsccStates, result); - - //now select the minimal/maximal element - if (minimize){ - result[state] = *std::min_element(reachableResults.begin(), reachableResults.end()); - } else { - result[state] = *std::max_element(reachableResults.begin(), reachableResults.end()); + storm::storage::BitVector statesInThisBsccs(numOfStates); + for (auto const& state : bscc) { + statesInThisBsccs.set(state); } + + ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); + + //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 reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBsccs, false); + + storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); + storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); } return result; } template - std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { - STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); - + std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType = boost::optional()) { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative))); } diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index d1e880e1a..d1f3edb60 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -37,7 +37,7 @@ namespace storm { std::vector computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const; std::vector computeCumulativeRewardsHelper(uint_fast64_t stepBound) const; std::vector computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const; - std::vector computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const; + std::vector computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const; static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc); From 9a83dfac10d24512640496a5cc65ce5d0ec44680 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 10 Apr 2015 18:33:26 +0200 Subject: [PATCH 07/20] Typo in DTMC, tried to use same approach for MDPs, which won't work. Former-commit-id: 5c1e835d09827fae0732b3b2eb4b1fc4f9b30cc6 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 6 +- .../prctl/SparseMdpPrctlModelChecker.cpp | 68 +++++++------------ 2 files changed, 29 insertions(+), 45 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 1f8642b54..56c4aeb5c 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -338,9 +338,9 @@ namespace storm { for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - storm::storage::BitVector statesInThisBsccs(numOfStates); + storm::storage::BitVector statesInThisBscc(numOfStates); for (auto const& state : bscc) { - statesInThisBsccs.set(state); + statesInThisBscc.set(state); } ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); @@ -349,7 +349,7 @@ namespace storm { //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 reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBsccs, false); + std::vector reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false); storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index cccc2bbe2..e27a179fb 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -339,61 +339,45 @@ namespace storm { // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); - // Now start with compute the long-run average for all end components in isolation. - std::vector lraValuesForEndComponents; - - // While doing so, we already gather some information for the following steps. - std::vector stateToMecIndexMap(numOfStates); - storm::storage::BitVector statesInMecs(numOfStates); + // First we check which states are in MECs. We use this later to speed up reachability analysis + storm::storage::BitVector statesNotInMecs(numOfStates, true); for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; // Gather information for later use. for (auto const& stateChoicesPair : mec) { - uint_fast64_t state = stateChoicesPair.first; - - statesInMecs.set(state); - stateToMecIndexMap[state] = currentMecIndex; + statesInMecs.set(stateChoicesPair.first, false); } - - // Compute the LRA value for the current MEC. - lraValuesForEndComponents.push_back(this->computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec)); } - // For fast transition rewriting, we build some auxiliary data structures. - storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; - // Prepare result vector. std::vector result(numOfStates); - //Set the values for all states in MECs. - for (auto state : statesInMecs) { - result[state] = lraValuesForEndComponents[stateToMecIndexMap[state]]; - } + //FIXME: This doesn't work for MDP. Simple counterexample, an MPD with precisely two MECs, both with LRA of 1. Max Reach Prob for both is 1 from the initial state + // Then this approach would yield an LRA value of 2 for the initial state, which is incorrect. - //for all states not in any mec set the result to the minimal/maximal value of the reachable MECs - //there might be a more efficient way to do this... - for (auto state : statesNotContainedInAnyMec){ - - //calculate what result values the reachable states in MECs have - storm::storage::BitVector currentState(numOfStates); - currentState.set(state); - storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( - transitionMatrix, currentState, storm::storage::BitVector(numOfStates, true), statesInMecs - ); - - storm::storage::BitVector reachableMecStates = statesInMecs & reachableStates; - std::vector reachableResults(reachableMecStates.getNumberOfSetBits()); - storm::utility::vector::selectVectorValues(reachableResults, reachableMecStates, result); - - //now select the minimal/maximal element - if (minimize){ - result[state] = *std::min_element(reachableResults.begin(), reachableResults.end()); - } else { - result[state] = *std::max_element(reachableResults.begin(), reachableResults.end()); - } - } + ////now we calculate the long run average for each MEC in isolation and compute the weighted contribution of the MEC to the LRA value of all states + //for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { + // storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; + + // storm::storage::BitVector statesInThisMec(numOfStates); + // for (auto const& state : bscc) { + // statesInThisMec.set(state); + // } + + // // Compute the LRA value for the current MEC. + // lraValuesForEndComponents.push_back(this->computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec)); + + // //the LRA value of a MEC contributes to the LRA value of a state with the probability of reaching that MEC from that state + // //thus we add Min/MaxProb(Eventually statesInThisBscc) * lraForThisBscc to the result vector + + // //the reachability probabilities will be zero in other MECs, thus we can set the left operand of the until formula to statesNotInMecs as an optimization + // std::vector reachProbThisBscc = this->computeUntilProbabilitiesHelper(minimize, statesNotInMecs, statesInThisMec, false); + + // storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); + // storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); + //} return result; } From 916c821b3e08c0fa2c0c8509aa21ec728e09e305 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 17 Apr 2015 16:55:12 +0200 Subject: [PATCH 08/20] Compute steady state for all BSCCs together by solving just one equation system instead of solving an equation system for each BSCC. Former-commit-id: 74f715c3a8c2c6e7dce03fb9acf736d517ed8bd1 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 40 ++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 56c4aeb5c..fc41e4d05 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -321,6 +321,7 @@ namespace storm { // First we check which states are in BSCCs. We use this later to speed up reachability analysis storm::storage::BitVector statesNotInBsccs(numOfStates, true); + std::vector stateToBsccIndexMap(model.getNumberOfStates()); for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; @@ -328,9 +329,45 @@ namespace storm { // Gather information for later use. for (auto const& state : bscc) { statesNotInBsccs.set(state, false); + stateToBsccIndexMap[state] = currentBsccIndex; } } + uint_fast64_t numStatesNotInBsccs = statesNotInBsccs.size() - statesNotInBsccs.getNumberOfSetBits(); + + // calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of the transposed transition matrix for the bsccs + std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); + storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, !statesNotInBsccs, !statesNotInBsccs, true).transpose(); + + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); + + for (index_type row = 0; row < this->getRowCount(); ++row) { + for (auto& entry : this->getRow(row)) { + if (entry.getColumn() == row) { + entry.setValue(entry.getValue() - one); + } + } + } + + std::vector bsccEquationSystemRightSide(numStatesNotInBsccs, zero); + std::vector bsccEquationSystemSolution(numStatesNotInBsccs, one); + solver->solveEquationSystem(bsccEquationSystem, 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 + vector bsccLra(bsccDecomposition.size(), zero); + vector bsccTotalValue(bsccDecomposition.size(), zero); + for (size_t i = 0, auto stateIter = !statesNotInBsccs.begin(); stateIter != !statesNotInBsccs.end(); ++i, ++stateIter) { + if (psiStates.get(*stateIter)) { + bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; + } + bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; + } + for (size_t i = 0; i < bsccLra.size(); ++i) { + bsccLra[i] /= bsccTotalValue[i]; + } + // Prepare result vector. std::vector result(numOfStates); @@ -343,7 +380,8 @@ namespace storm { statesInThisBscc.set(state); } - ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); + //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 From 0fdb3685d1f9fafe518dded4f4399cab28ddad4a Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 17 Apr 2015 18:02:41 +0200 Subject: [PATCH 09/20] Computing LRA for states not in bsccs as expected reward Former-commit-id: 4bcb5f0a6e0044b40b0d74fe820a73ebee691636 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 104 +++++++++++++----- 1 file changed, 74 insertions(+), 30 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index fc41e4d05..2ef77454f 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -320,80 +320,124 @@ namespace storm { typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); // First we check which states are in BSCCs. We use this later to speed up reachability analysis - storm::storage::BitVector statesNotInBsccs(numOfStates, true); - std::vector stateToBsccIndexMap(model.getNumberOfStates()); + storm::storage::BitVector statesInBsccs(numOfStates); + + std::vector 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) { - statesNotInBsccs.set(state, false); + statesInBsccs.set(state); stateToBsccIndexMap[state] = currentBsccIndex; } } - uint_fast64_t numStatesNotInBsccs = statesNotInBsccs.size() - statesNotInBsccs.getNumberOfSetBits(); + 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 std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); - storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, !statesNotInBsccs, !statesNotInBsccs, true).transpose(); + storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true).transpose(); ValueType one = storm::utility::one(); ValueType zero = storm::utility::zero(); - for (index_type row = 0; row < this->getRowCount(); ++row) { - for (auto& entry : this->getRow(row)) { + 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); } } } - std::vector bsccEquationSystemRightSide(numStatesNotInBsccs, zero); - std::vector bsccEquationSystemSolution(numStatesNotInBsccs, one); + std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); + std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); solver->solveEquationSystem(bsccEquationSystem, 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 - vector bsccLra(bsccDecomposition.size(), zero); - vector bsccTotalValue(bsccDecomposition.size(), zero); - for (size_t i = 0, auto stateIter = !statesNotInBsccs.begin(); stateIter != !statesNotInBsccs.end(); ++i, ++stateIter) { + std::vector bsccLra(bsccDecomposition.size(), zero); + std::vector 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 (size_t i = 0; i < bsccLra.size(); ++i) { + for (i = 0; i < bsccLra.size(); ++i) { bsccLra[i] /= bsccTotalValue[i]; } - // Prepare result vector. - std::vector result(numOfStates); + //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 - //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]; + std::vector rewardRightSide; + rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); - storm::storage::BitVector statesInThisBscc(numOfStates); - for (auto const& state : bscc) { - statesInThisBscc.set(state); + 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 rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs); + rewardEquationSystemMatrix.convertToEquationSystem(); + + std::vector rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one); - //ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); - ValueType lraForThisBscc = bsccLra[currentBsccIndex]; + solver->solveEquationSystem(rewardEquationSystemMatrix, rewardSolution, rewardRightSide); - //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 + // now fill the result vector + std::vector result(numOfStates); - //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 reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false); - - storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); - storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); + 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 reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false); + // + // storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); + // storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); + //} + + //return result; } template From e4968b1ddee0819fd399366e321398f39dcc66ea Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 24 Apr 2015 08:59:34 +0200 Subject: [PATCH 10/20] Fixed minor issue in cli Former-commit-id: ed63925765551ec91d8db0262e6797c4de6d3fe8 --- examples/dtmc/crowds/crowds15_5.pm | 6 +++--- examples/dtmc/crowds/crowds20_5.pm | 7 ++++--- src/utility/cli.h | 2 +- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/dtmc/crowds/crowds15_5.pm b/examples/dtmc/crowds/crowds15_5.pm index 511d962f0..95ab8a910 100644 --- a/examples/dtmc/crowds/crowds15_5.pm +++ b/examples/dtmc/crowds/crowds15_5.pm @@ -2,9 +2,9 @@ dtmc // probability of forwarding const double PF = 0.8; -const double notPF = .2; // must be 1-PF +const double notPF = 0.2; // must be 1-PF // probability that a crowd member is bad -const double badC = .167; +const double badC = 0.167; // probability that a crowd member is good const double goodC = 0.833; // Total number of protocol runs to analyze @@ -92,4 +92,4 @@ endmodule label "observe0Greater1" = observe0 > 1; label "observeIGreater1" = observe1 > 1 | observe2 > 1 | observe3 > 1 | observe4 > 1 | observe5 > 1 | observe6 > 1 | observe7 > 1 | observe8 > 1 | observe9 > 1 | observe10 > 1 | observe11 > 1 | observe12 > 1 | observe13 > 1 | observe14 > 1; -label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1; \ No newline at end of file +label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1; diff --git a/examples/dtmc/crowds/crowds20_5.pm b/examples/dtmc/crowds/crowds20_5.pm index 31c63770a..1809c22d7 100644 --- a/examples/dtmc/crowds/crowds20_5.pm +++ b/examples/dtmc/crowds/crowds20_5.pm @@ -2,9 +2,9 @@ dtmc // probability of forwarding const double PF = 0.8; -const double notPF = .2; // must be 1-PF +const double notPF = 0.2; // must be 1-PF // probability that a crowd member is bad -const double badC = .167; +const double badC = 0.167; // probability that a crowd member is good const double goodC = 0.833; // Total number of protocol runs to analyze @@ -107,4 +107,5 @@ endmodule label "observe0Greater1" = observe0 > 1; label "observeIGreater1" = observe1 > 1 | observe2 > 1 | observe3 > 1 | observe4 > 1 | observe5 > 1 | observe6 > 1 | observe7 > 1 | observe8 > 1 | observe9 > 1 | observe10 > 1 | observe11 > 1 | observe12 > 1 | observe13 > 1 | observe14 > 1 | observe15 > 1 | observe16 > 1 | observe17 > 1 | observe18 > 1 | observe19 > 1; -label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1 & observe15 <= 1 & observe16 <= 1 & observe17 <= 1 & observe18 <= 1 & observe19 <= 1; \ No newline at end of file +label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1 & observe15 <= 1 & observe16 <= 1 & observe17 <= 1 & observe18 <= 1 & observe19 <= 1; + diff --git a/src/utility/cli.h b/src/utility/cli.h index 780201cf3..7fd03536d 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -426,7 +426,7 @@ namespace storm { } else { storm::modelchecker::SparseDtmcEliminationModelChecker modelchecker2(*dtmc); if (modelchecker2.canHandle(*formula.get())) { - modelchecker2.check(*formula.get()); + result = modelchecker2.check(*formula.get()); } } } else if (model->getType() == storm::models::ModelType::Mdp) { From 8fc58439bcc59a66b7f1b4e4ad0c10081f12f331 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Tue, 28 Apr 2015 19:56:02 +0200 Subject: [PATCH 11/20] Computing LRA as expected reward in MDPs. - Everything compiles without error. No tests yet. Former-commit-id: d8cceb02fcffa27fb91acf6b5136fc33fb45865c --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 2 +- .../prctl/SparseMdpPrctlModelChecker.cpp | 81 ++++++++++++------- 2 files changed, 54 insertions(+), 29 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 2ef77454f..2661a8d74 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -314,7 +314,7 @@ namespace storm { } // Start by decomposing the DTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(this->getModelAs>(), false, true); + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(this->getModel(), false, true); // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index e27a179fb..e45cc1baa 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -334,50 +334,75 @@ namespace storm { } // Start by decomposing the MDP into its MECs. - storm::storage::MaximalEndComponentDecomposition mecDecomposition(this->getModelAs>()); + storm::storage::MaximalEndComponentDecomposition mecDecomposition(this->getModel()); // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); - // First we check which states are in MECs. We use this later to speed up reachability analysis - storm::storage::BitVector statesNotInMecs(numOfStates, true); + //first calculate LRA for the Maximal End Components. + storm::storage::BitVector statesInMecs(numOfStates); + std::vector stateToMecIndexMap(transitionMatrix.getColumnCount()); + std::vector mecLra(mecDecomposition.size(), zero); for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; + mecLra[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec); + // Gather information for later use. for (auto const& stateChoicesPair : mec) { - statesInMecs.set(stateChoicesPair.first, false); + statesInMecs.set(stateChoicesPair.first); + stateToMecIndexMap[stateChoicesPair.first] = currentMecIndex; } } - // Prepare result vector. - std::vector result(numOfStates); - - //FIXME: This doesn't work for MDP. Simple counterexample, an MPD with precisely two MECs, both with LRA of 1. Max Reach Prob for both is 1 from the initial state - // Then this approach would yield an LRA value of 2 for the initial state, which is incorrect. - - ////now we calculate the long run average for each MEC in isolation and compute the weighted contribution of the MEC to the LRA value of all states - //for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { - // storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; - - // storm::storage::BitVector statesInThisMec(numOfStates); - // for (auto const& state : bscc) { - // statesInThisMec.set(state); - // } - - // // Compute the LRA value for the current MEC. - // lraValuesForEndComponents.push_back(this->computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec)); + //calculate LRA for states not in MECs as expected reachability rewards + //we add an auxiliary target state, every state in any MEC has a choice to move to that state with prob 1. + //This transitions have the LRA of the MEC as reward. + //The expected reward corresponds to sum of LRAs in MEC weighted by the reachability probability of the MEC + + std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); + + //we now build the submatrix of the transition matrix of the system with the auxiliary state, that only contains the states from + //the original state, i.e. all "maybe-states" + storm::storage::SparseMatrixBuilder rewardEquationSystemMatrixBuilder(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), + transitionMatrix.getColumnCount(), + transitionMatrix.getEntryCount(), + true, + true, + transitionMatrix.getRowGroupCount()); + + std::vector rewardRightSide(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), zero); + + uint_fast64_t rowIndex = 0; + uint_fast64_t oldRowIndex = 0; + for (uint_fast64_t rowGroupIndex = 0; rowGroupIndex < transitionMatrix.getRowGroupCount(); ++rowGroupIndex) { + rewardEquationSystemMatrixBuilder.newRowGroup(rowIndex); + for (uint_fast64_t i = 0; i < transitionMatrix.getRowGroupSize(rowGroupIndex); ++i) { + for (auto entry : transitionMatrix.getRow(oldRowIndex)) { + //copy over values from transition matrix of the actual system + rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, entry.getColumn(), entry.getValue()); + } + ++oldRowIndex; + ++rowIndex; + } + ++rowGroupIndex; + if (statesInMecs.get(rowGroupIndex)) { + //add the choice where we go to the auxiliary state, which is a row with all zeros in the submatrix we build + ++rowIndex; + //put the transition-reward on the right side + rewardRightSide[rowIndex] = mecLra[rowGroupIndex]; + } + } - // //the LRA value of a MEC contributes to the LRA value of a state with the probability of reaching that MEC from that state - // //thus we add Min/MaxProb(Eventually statesInThisBscc) * lraForThisBscc to the result vector + storm::storage::SparseMatrix rewardEquationSystemMatrix = rewardEquationSystemMatrixBuilder.build(); + rewardEquationSystemMatrix.convertToEquationSystem(); - // //the reachability probabilities will be zero in other MECs, thus we can set the left operand of the until formula to statesNotInMecs as an optimization - // std::vector reachProbThisBscc = this->computeUntilProbabilitiesHelper(minimize, statesNotInMecs, statesInThisMec, false); + std::vector result(rewardEquationSystemMatrix.getColumnCount(), one); - // storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); - // storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); - //} + solver->solveEquationSystem(rewardEquationSystemMatrix, result, rewardRightSide); return result; } From 1f87e7c8b269027c7a4de105bcd9324933d90d66 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Wed, 29 Apr 2015 20:35:16 +0200 Subject: [PATCH 12/20] First test for LRA on MDPs Former-commit-id: c022ceaf01964fa47ae4e85173673714e2879ec9 --- .../SparseMdpPrctlModelCheckerTest.cpp | 99 ++++++++++++------- 1 file changed, 64 insertions(+), 35 deletions(-) diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index e7fb71382..b6483b306 100644 --- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp @@ -143,54 +143,83 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) { ASSERT_EQ(3172ull, mdp->getNumberOfStates()); ASSERT_EQ(7144ull, mdp->getNumberOfTransitions()); - storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::shared_ptr>(new storm::solver::NativeNondeterministicLinearEquationSolver())); + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::shared_ptr>(new storm::solver::NativeNondeterministicLinearEquationSolver())); - auto labelFormula = std::make_shared("elected"); - auto eventuallyFormula = std::make_shared(labelFormula); - auto minProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, eventuallyFormula); + auto labelFormula = std::make_shared("elected"); + auto eventuallyFormula = std::make_shared(labelFormula); + auto minProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, eventuallyFormula); + + std::unique_ptr result = checker.check(*minProbabilityOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); - std::unique_ptr result = checker.check(*minProbabilityOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); - EXPECT_NEAR(1, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - auto maxProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, eventuallyFormula); - - result = checker.check(*maxProbabilityOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); - + auto maxProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, eventuallyFormula); + + result = checker.check(*maxProbabilityOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(1, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - labelFormula = std::make_shared("elected"); - auto trueFormula = std::make_shared(true); - auto boundedUntilFormula = std::make_shared(trueFormula, labelFormula, 25); - minProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, boundedUntilFormula); - - result = checker.check(*minProbabilityOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); - + labelFormula = std::make_shared("elected"); + auto trueFormula = std::make_shared(true); + auto boundedUntilFormula = std::make_shared(trueFormula, labelFormula, 25); + minProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, boundedUntilFormula); + + result = checker.check(*minProbabilityOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.0625, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - maxProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, boundedUntilFormula); - - result = checker.check(*maxProbabilityOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); - + maxProbabilityOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, boundedUntilFormula); + + result = checker.check(*maxProbabilityOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.0625, quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - labelFormula = std::make_shared("elected"); - auto reachabilityRewardFormula = std::make_shared(labelFormula); - auto minRewardOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula); + labelFormula = std::make_shared("elected"); + auto reachabilityRewardFormula = std::make_shared(labelFormula); + auto minRewardOperatorFormula = std::make_shared(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula); + + result = checker.check(*minRewardOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult5 = result->asExplicitQuantitativeCheckResult(); - result = checker.check(*minRewardOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult5 = result->asExplicitQuantitativeCheckResult(); - EXPECT_NEAR(4.285689611, quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - auto maxRewardOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula); + auto maxRewardOperatorFormula = std::make_shared(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula); + + result = checker.check(*maxRewardOperatorFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult6 = result->asExplicitQuantitativeCheckResult(); - result = checker.check(*maxRewardOperatorFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult6 = result->asExplicitQuantitativeCheckResult(); - EXPECT_NEAR(4.285689611, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision()); } + +TEST(SparseMdpPrctlModelCheckerTest, LRA) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> mdp; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(2, 2, 2); + matrixBuilder.addNextValue(0, 1, 1.); + matrixBuilder.addNextValue(1, 0, 1.); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::AtomicPropositionsLabeling ap(2, 1); + ap.addAtomicProposition("a"); + ap.addAtomicPropositionToState("a", 1); + + mdp.reset(new storm::models::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::shared_ptr>(new storm::solver::NativeNondeterministicLinearEquationSolver())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + std::unique_ptr result = checker.check(*lraFormula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } +} From d4f051c4f0e5e799c9e9090c316f5b7a6f719669 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Fri, 1 May 2015 14:46:26 +0200 Subject: [PATCH 13/20] Fixed Windows build Former-commit-id: 53c99736de67fcf7020ec043c2c03dbe87e669b9 --- .../cudd-2.5.0/src/cudd/cuddAddApply.c | 6 +++--- .../GmmxxCtmcCslModelCheckerTest.cpp | 21 +++++++++++++++---- .../GmmxxHybridDtmcPrctlModelCheckerTest.cpp | 12 +++++++++-- .../SparseCtmcCslModelCheckerTest.cpp | 18 +++++++++++++--- 4 files changed, 45 insertions(+), 12 deletions(-) diff --git a/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c b/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c index fc1ab8c15..e789bfc73 100644 --- a/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c +++ b/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c @@ -1078,7 +1078,7 @@ Cudd_addMod( F = *f; G = *g; if (cuddIsConstant(F) && cuddIsConstant(G)) { // If g is <=0, then result is NaN - if (cuddV(G) <= 0) value = (0.0/0.0); + if (cuddV(G) <= 0) value = (NAN); // Take care of negative case (% is remainder, not modulo) else { rem = ((int)cuddV(F) % (int)cuddV(G)); @@ -1119,9 +1119,9 @@ Cudd_addLogXY( F = *f; G = *g; if (cuddIsConstant(F) && cuddIsConstant(G)) { // If base is <=0 or ==1 (or +Inf/NaN), then result is NaN - if (cuddV(G) <= 0 || cuddV(G) == 1.0 || G==DD_PLUS_INFINITY(dd) || cuddV(G) != cuddV(G)) value = (0.0/0.0); + if (cuddV(G) <= 0 || cuddV(G) == 1.0 || G==DD_PLUS_INFINITY(dd) || cuddV(G) != cuddV(G)) value = (NAN); // If arg is <0 or NaN, then result is NaN - else if (cuddV(F) < 0 || cuddV(F) != cuddV(F)) value = (0.0/0.0); + else if (cuddV(F) < 0 || cuddV(F) != cuddV(F)) value = (NAN); // If arg is +Inf, then result is +Inf else if (F==DD_PLUS_INFINITY(dd)) return DD_PLUS_INFINITY(dd); // If arg is (positive/negative) 0, then result is -Inf diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp index 6a814d63d..5c637afcb 100644 --- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp @@ -22,7 +22,11 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) { std::shared_ptr formula(nullptr); // Build the model. - typename storm::builder::ExplicitPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "num_repairs"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); @@ -80,7 +84,12 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) { std::shared_ptr formula(nullptr); // Build the model. - typename storm::builder::ExplicitPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif + options.buildRewards = true; options.rewardModelName = "up"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); @@ -172,8 +181,12 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) { std::shared_ptr formula(nullptr); // Build the model. - typename storm::builder::ExplicitPrismModelBuilder::Options options; - options.buildRewards = true; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif + options.buildRewards = true; options.rewardModelName = "customers"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType()); diff --git a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp index cd1a014a9..df17ce336 100644 --- a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp @@ -16,7 +16,11 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, 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. - typename storm::builder::DdPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::DdPrismModelBuilder::Options options; +#else + typename storm::builder::DdPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "coin_flips"; std::shared_ptr> model = storm::builder::DdPrismModelBuilder::translateProgram(program, options); @@ -118,7 +122,11 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, 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. - typename storm::builder::DdPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::DdPrismModelBuilder::Options options; +#else + typename storm::builder::DdPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "num_rounds"; std::shared_ptr> model = storm::builder::DdPrismModelBuilder::translateProgram(program, options); diff --git a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp index 3dfa1b6be..9db619115 100644 --- a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp @@ -22,7 +22,11 @@ TEST(SparseCtmcCslModelCheckerTest, Cluster) { std::shared_ptr formula(nullptr); // Build the model. - typename storm::builder::ExplicitPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "num_repairs"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); @@ -80,7 +84,11 @@ TEST(SparseCtmcCslModelCheckerTest, Embedded) { std::shared_ptr formula(nullptr); // Build the model. - typename storm::builder::ExplicitPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "up"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); @@ -172,7 +180,11 @@ TEST(SparseCtmcCslModelCheckerTest, Tandem) { std::shared_ptr formula(nullptr); // Build the model with the customers reward structure. - typename storm::builder::ExplicitPrismModelBuilder::Options options; +#ifdef WINDOWS + storm::builder::ExplicitPrismModelBuilder::Options options; +#else + typename storm::builder::ExplicitPrismModelBuilder::Options options; +#endif options.buildRewards = true; options.rewardModelName = "customers"; std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); From 716cf3abdd2c0fc99a8f1cc5b017ae22d5f526b7 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Mon, 4 May 2015 23:35:35 +0200 Subject: [PATCH 14/20] Adapted to new solver interface some tests and bugfixes. Tests still failing. Former-commit-id: da3b75aefd5d5a4c6f158f0c4aecb9e58dd7fc39 --- src/modelchecker/AbstractModelChecker.cpp | 8 ++--- .../prctl/SparseDtmcPrctlModelChecker.cpp | 20 +++++++++---- .../prctl/SparseDtmcPrctlModelChecker.h | 2 +- .../prctl/SparseMdpPrctlModelChecker.cpp | 29 ++++++++++++++----- .../SparseMdpPrctlModelCheckerTest.cpp | 14 ++++----- 5 files changed, 48 insertions(+), 25 deletions(-) diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp index 7ba1de099..2b6d2ef3a 100644 --- a/src/modelchecker/AbstractModelChecker.cpp +++ b/src/modelchecker/AbstractModelChecker.cpp @@ -216,8 +216,8 @@ namespace storm { } } - std::unique_ptr AbstractModelChecker::checkExpectedTimeOperatorFormula(storm::logic::ExpectedTimeOperatorFormula const& stateFormula) { - STORM_LOG_THROW(stateFormula.getSubformula().isStateFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid."); + std::unique_ptr AbstractModelChecker::checkExpectedTimeOperatorFormula(storm::logic::ExpectedTimeOperatorFormula const& stateFormula) { + STORM_LOG_THROW(stateFormula.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid."); // If the reward bound is 0, is suffices to do qualitative model checking. bool qualitative = false; @@ -248,8 +248,8 @@ namespace storm { } } - std::unique_ptr AbstractModelChecker::checkLongRunAverageOperatorFormula(storm::logic::LongRunAverageOperatorFormula const& stateFormula) { - STORM_LOG_THROW(stateFormula.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid."); + std::unique_ptr AbstractModelChecker::checkLongRunAverageOperatorFormula(storm::logic::LongRunAverageOperatorFormula const& stateFormula) { + STORM_LOG_THROW(stateFormula.getSubformula().isStateFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid."); std::unique_ptr result; if (stateFormula.hasOptimalityType()) { diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 02a2b5758..9c70023a7 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -340,7 +340,6 @@ namespace storm { 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 - std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true).transpose(); ValueType one = storm::utility::one(); @@ -356,7 +355,10 @@ namespace storm { std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); - solver->solveEquationSystem(bsccEquationSystem, bsccEquationSystemSolution, bsccEquationSystemRightSide); + { + 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 @@ -395,7 +397,10 @@ namespace storm { std::vector rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one); - solver->solveEquationSystem(rewardEquationSystemMatrix, rewardSolution, rewardRightSide); + { + auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix); + solver->solveEquationSystem(rewardSolution, rewardRightSide); + } // now fill the result vector std::vector result(numOfStates); @@ -453,7 +458,7 @@ namespace storm { template - ValueType SparseDtmcPrctlModelChecker::computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc) { + ValueType SparseDtmcPrctlModelChecker::computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { //if size is one we can immediately derive the result if (bscc.size() == 1){ if (psiStates.get(*(bscc.begin()))) { @@ -462,7 +467,6 @@ namespace storm { return storm::utility::zero(); } } - std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount()); subsystem.set(bscc.begin(), bscc.end()); @@ -502,7 +506,11 @@ namespace storm { std::vector b(subsystemMatrix.getRowCount(), zero); b[subsystemMatrix.getRowCount() - 1] = one; - solver->solveEquationSystem(subsystemMatrix, steadyStateDist, b); + + { + 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(); diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index a84413495..b55f00faa 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -43,7 +43,7 @@ namespace storm { static std::vector computeReachabilityRewardsHelper(storm::storage::SparseMatrix const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative); std::vector computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const; - static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc); + static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); // An object that is used for retrieving linear equation solvers. std::unique_ptr> linearEquationSolverFactory; diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index ba88c4ea7..9fa78f041 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -366,14 +366,12 @@ namespace storm { //This transitions have the LRA of the MEC as reward. //The expected reward corresponds to sum of LRAs in MEC weighted by the reachability probability of the MEC - std::unique_ptr> solver = storm::utility::solver::getLinearEquationSolver(); - //we now build the submatrix of the transition matrix of the system with the auxiliary state, that only contains the states from //the original state, i.e. all "maybe-states" storm::storage::SparseMatrixBuilder rewardEquationSystemMatrixBuilder(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), transitionMatrix.getColumnCount(), transitionMatrix.getEntryCount(), - true, + false, true, transitionMatrix.getRowGroupCount()); @@ -384,19 +382,32 @@ namespace storm { for (uint_fast64_t rowGroupIndex = 0; rowGroupIndex < transitionMatrix.getRowGroupCount(); ++rowGroupIndex) { rewardEquationSystemMatrixBuilder.newRowGroup(rowIndex); for (uint_fast64_t i = 0; i < transitionMatrix.getRowGroupSize(rowGroupIndex); ++i) { + //we have to make sure that an entry exists for all diagonal elements, even if it is zero. Other wise the call to convertToEquationSystem will produce wrong results or fail. + bool foundDiagonal = false; for (auto entry : transitionMatrix.getRow(oldRowIndex)) { + if (!foundDiagonal) { + if (entry.getColumn() > rowGroupIndex) { + foundDiagonal = true; + rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); + } else if (entry.getColumn() == rowGroupIndex) { + foundDiagonal = true; + } + } //copy over values from transition matrix of the actual system rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, entry.getColumn(), entry.getValue()); } + if (!foundDiagonal) { + rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); + } ++oldRowIndex; ++rowIndex; } - ++rowGroupIndex; if (statesInMecs.get(rowGroupIndex)) { + //put the transition-reward on the right side + rewardRightSide[rowIndex] = mecLra[stateToMecIndexMap[rowGroupIndex]]; //add the choice where we go to the auxiliary state, which is a row with all zeros in the submatrix we build + rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); ++rowIndex; - //put the transition-reward on the right side - rewardRightSide[rowIndex] = mecLra[rowGroupIndex]; } } @@ -405,7 +416,11 @@ namespace storm { std::vector result(rewardEquationSystemMatrix.getColumnCount(), one); - solver->solveEquationSystem(rewardEquationSystemMatrix, result, rewardRightSide); + { + auto solver = this->MinMaxLinearEquationSolverFactory->create(rewardEquationSystemMatrix); + solver->solveEquationSystem(minimize, result, rewardRightSide); + } + return result; } diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index 9a1ab6b1d..3ea26261f 100644 --- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp @@ -197,7 +197,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) { TEST(SparseMdpPrctlModelCheckerTest, LRA) { storm::storage::SparseMatrixBuilder matrixBuilder; - std::shared_ptr> mdp; + std::shared_ptr> mdp; { matrixBuilder = storm::storage::SparseMatrixBuilder(2, 2, 2); @@ -205,16 +205,16 @@ TEST(SparseMdpPrctlModelCheckerTest, LRA) { matrixBuilder.addNextValue(1, 0, 1.); storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); - storm::models::AtomicPropositionsLabeling ap(2, 1); - ap.addAtomicProposition("a"); - ap.addAtomicPropositionToState("a", 1); + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); - mdp.reset(new storm::models::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + mdp.reset(new storm::models::sparse::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::shared_ptr>(new storm::solver::NativeNondeterministicLinearEquationSolver())); + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory())); auto labelFormula = std::make_shared("a"); - auto lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); std::unique_ptr result = checker.check(*lraFormula); storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult = result->asExplicitQuantitativeCheckResult(); From dfab1c291ce464736ee6f3fc0f1fe1a4f69c20a1 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Wed, 6 May 2015 23:05:38 +0200 Subject: [PATCH 15/20] Error fixed. Former-commit-id: b2514633a75876f8fcca44d53a1a35252f78d4f0 --- src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp index abec6a929..ebedb6c25 100644 --- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp @@ -390,7 +390,7 @@ namespace storm { // Set the values for all states in MECs. for (auto state : statesInMecs) { - result[state] = lraValuesForEndComponents[stateToMecIndexMap[state]]; + result[state] = x[firstAuxiliaryStateIndex + stateToMecIndexMap[state]]; } return result; From 0ba629ad3fa0979678cf8bb9c21d174bed291eef Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Wed, 6 May 2015 23:06:01 +0200 Subject: [PATCH 16/20] More tests, bugfixes: All tests pass. Former-commit-id: f37c02a9d773c91e7a3dc72c3dbc19a907242661 --- .../prctl/SparseMdpPrctlModelChecker.cpp | 169 ++++++---- .../SparseMdpPrctlModelCheckerTest.cpp | 314 +++++++++++++++++- 2 files changed, 422 insertions(+), 61 deletions(-) diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 9fa78f041..00045105d 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -341,18 +341,19 @@ namespace storm { // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); + std::vector const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); ValueType one = storm::utility::one(); ValueType zero = storm::utility::zero(); //first calculate LRA for the Maximal End Components. storm::storage::BitVector statesInMecs(numOfStates); std::vector stateToMecIndexMap(transitionMatrix.getColumnCount()); - std::vector mecLra(mecDecomposition.size(), zero); + std::vector lraValuesForEndComponents(mecDecomposition.size(), zero); for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; - mecLra[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec); + lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec); // Gather information for later use. for (auto const& stateChoicesPair : mec) { @@ -361,66 +362,122 @@ namespace storm { } } - //calculate LRA for states not in MECs as expected reachability rewards - //we add an auxiliary target state, every state in any MEC has a choice to move to that state with prob 1. - //This transitions have the LRA of the MEC as reward. - //The expected reward corresponds to sum of LRAs in MEC weighted by the reachability probability of the MEC - - //we now build the submatrix of the transition matrix of the system with the auxiliary state, that only contains the states from - //the original state, i.e. all "maybe-states" - storm::storage::SparseMatrixBuilder rewardEquationSystemMatrixBuilder(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), - transitionMatrix.getColumnCount(), - transitionMatrix.getEntryCount(), - false, - true, - transitionMatrix.getRowGroupCount()); - - std::vector rewardRightSide(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), zero); - - uint_fast64_t rowIndex = 0; - uint_fast64_t oldRowIndex = 0; - for (uint_fast64_t rowGroupIndex = 0; rowGroupIndex < transitionMatrix.getRowGroupCount(); ++rowGroupIndex) { - rewardEquationSystemMatrixBuilder.newRowGroup(rowIndex); - for (uint_fast64_t i = 0; i < transitionMatrix.getRowGroupSize(rowGroupIndex); ++i) { - //we have to make sure that an entry exists for all diagonal elements, even if it is zero. Other wise the call to convertToEquationSystem will produce wrong results or fail. - bool foundDiagonal = false; - for (auto entry : transitionMatrix.getRow(oldRowIndex)) { - if (!foundDiagonal) { - if (entry.getColumn() > rowGroupIndex) { - foundDiagonal = true; - rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); - } else if (entry.getColumn() == rowGroupIndex) { - foundDiagonal = true; - } + // For fast transition rewriting, we build some auxiliary data structures. + storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; + uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits(); + uint_fast64_t lastStateNotInMecs = 0; + uint_fast64_t numberOfStatesNotInMecs = 0; + std::vector statesNotInMecsBeforeIndex; + statesNotInMecsBeforeIndex.reserve(this->getModel().getNumberOfStates()); + for (auto state : statesNotContainedInAnyMec) { + while (lastStateNotInMecs <= state) { + statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs); + ++lastStateNotInMecs; + } + ++numberOfStatesNotInMecs; + } + + // Finally, we are ready to create the SSP matrix and right-hand side of the SSP. + std::vector b; + typename storm::storage::SparseMatrixBuilder sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size()); + + // If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications). + uint_fast64_t currentChoice = 0; + for (auto state : statesNotContainedInAnyMec) { + sspMatrixBuilder.newRowGroup(currentChoice); + + for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) { + std::vector auxiliaryStateToProbabilityMap(mecDecomposition.size()); + b.push_back(storm::utility::zero()); + + for (auto element : transitionMatrix.getRow(choice)) { + if (statesNotContainedInAnyMec.get(element.getColumn())) { + // If the target state is not contained in an MEC, we can copy over the entry. + sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue()); + } else { + // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector + // so that we are able to write the cumulative probability to the MEC into the matrix. + auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue(); } - //copy over values from transition matrix of the actual system - rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, entry.getColumn(), entry.getValue()); } - if (!foundDiagonal) { - rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); + + // Now insert all (cumulative) probability values that target an MEC. + for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) { + if (auxiliaryStateToProbabilityMap[mecIndex] != 0) { + sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]); + } } - ++oldRowIndex; - ++rowIndex; } - if (statesInMecs.get(rowGroupIndex)) { - //put the transition-reward on the right side - rewardRightSide[rowIndex] = mecLra[stateToMecIndexMap[rowGroupIndex]]; - //add the choice where we go to the auxiliary state, which is a row with all zeros in the submatrix we build - rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero); - ++rowIndex; + } + + // Now we are ready to construct the choices for the auxiliary states. + for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) { + storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex]; + sspMatrixBuilder.newRowGroup(currentChoice); + + for (auto const& stateChoicesPair : mec) { + uint_fast64_t state = stateChoicesPair.first; + boost::container::flat_set const& choicesInMec = stateChoicesPair.second; + + for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) { + std::vector auxiliaryStateToProbabilityMap(mecDecomposition.size()); + + // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state. + if (choicesInMec.find(choice) == choicesInMec.end()) { + b.push_back(storm::utility::zero()); + + for (auto element : transitionMatrix.getRow(choice)) { + if (statesNotContainedInAnyMec.get(element.getColumn())) { + // If the target state is not contained in an MEC, we can copy over the entry. + sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue()); + } else { + // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector + // so that we are able to write the cumulative probability to the MEC into the matrix. + auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue(); + } + } + + // Now insert all (cumulative) probability values that target an MEC. + for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) { + if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) { + // If the target MEC is the same as the current one, instead of adding a transition, we need to add the weighted reward + // to the right-hand side vector of the SSP. + if (mecIndex == targetMecIndex) { + b.back() += auxiliaryStateToProbabilityMap[mecIndex] * lraValuesForEndComponents[mecIndex]; + } else { + // Otherwise, we add a transition to the auxiliary state that is associated with the target MEC. + sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]); + } + } + } + + ++currentChoice; + } + } } + + // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC. + ++currentChoice; + b.push_back(lraValuesForEndComponents[mecIndex]); } - storm::storage::SparseMatrix rewardEquationSystemMatrix = rewardEquationSystemMatrixBuilder.build(); - rewardEquationSystemMatrix.convertToEquationSystem(); + // Finalize the matrix and solve the corresponding system of equations. + storm::storage::SparseMatrix sspMatrix = sspMatrixBuilder.build(currentChoice); - std::vector result(rewardEquationSystemMatrix.getColumnCount(), one); + std::vector sspResult(numberOfStatesNotInMecs + mecDecomposition.size()); + std::unique_ptr> solver = MinMaxLinearEquationSolverFactory->create(sspMatrix); + solver->solveEquationSystem(minimize, sspResult, b); - { - auto solver = this->MinMaxLinearEquationSolverFactory->create(rewardEquationSystemMatrix); - solver->solveEquationSystem(minimize, result, rewardRightSide); + // Prepare result vector. + std::vector result(this->getModel().getNumberOfStates()); + + // Set the values for states not contained in MECs. + storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, sspResult); + + // Set the values for all states in MECs. + for (auto state : statesInMecs) { + result[state] = sspResult[firstAuxiliaryStateIndex + stateToMecIndexMap[state]]; } - return result; } @@ -455,16 +512,16 @@ namespace storm { // Now, based on the type of the state, create a suitable constraint. for (auto choice : stateChoicesPair.second) { - storm::expressions::Expression constraint = solver->getConstant(1); - ValueType w = 0; + storm::expressions::Expression constraint = -lambda; + ValueType r = 0; for (auto element : transitionMatrix.getRow(choice)) { constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); if (psiStates.get(element.getColumn())) { - w += element.getValue(); + r += element.getValue(); } } - constraint = constraint - solver->getConstant(w) * lambda; + constraint = solver->getConstant(r) + constraint; if (minimize) { constraint = stateToVariableMap.at(state) <= constraint; diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index 3ea26261f..d3210a395 100644 --- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp @@ -195,7 +195,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) { EXPECT_NEAR(4.285689611, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision()); } -TEST(SparseMdpPrctlModelCheckerTest, LRA) { +TEST(SparseMdpPrctlModelCheckerTest, LRA_SingleMec) { storm::storage::SparseMatrixBuilder matrixBuilder; std::shared_ptr> mdp; @@ -216,10 +216,314 @@ TEST(SparseMdpPrctlModelCheckerTest, LRA) { auto labelFormula = std::make_shared("a"); auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); - std::unique_ptr result = checker.check(*lraFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult = result->asExplicitQuantitativeCheckResult(); + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + { + matrixBuilder = storm::storage::SparseMatrixBuilder(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 transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + mdp.reset(new storm::models::sparse::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(4, 3, 4, true, true, 3); + matrixBuilder.newRowGroup(0); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.newRowGroup(1); + matrixBuilder.addNextValue(1, 0, 1); + matrixBuilder.addNextValue(2, 2, 1); + matrixBuilder.newRowGroup(3); + matrixBuilder.addNextValue(3, 0, 1); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 2); + ap.addLabel("b"); + ap.addLabelToState("b", 0); + ap.addLabel("c"); + ap.addLabelToState("c", 0); + ap.addLabelToState("c", 2); + + mdp.reset(new storm::models::sparse::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + 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()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared("b"); + lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.5, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.5, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.5, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1. / 3., quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared("c"); + lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult5 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(2. / 3., quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(2. / 3., quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(2. / 3., quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(.5, quantitativeResult[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(.5, quantitativeResult[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult6 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.5, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.5, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.5, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision()); } } + +TEST(SparseMdpPrctlModelCheckerTest, LRA) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> mdp; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(4, 3, 4, true, true, 3); + matrixBuilder.newRowGroup(0); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.newRowGroup(1); + matrixBuilder.addNextValue(1, 1, 1); + matrixBuilder.addNextValue(2, 2, 1); + matrixBuilder.newRowGroup(3); + matrixBuilder.addNextValue(3, 2, 1); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 0); + ap.addLabel("b"); + ap.addLabelToState("b", 1); + ap.addLabel("c"); + ap.addLabelToState("c", 2); + + mdp.reset(new storm::models::sparse::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.0, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared("b"); + lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.0, quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + labelFormula = std::make_shared("c"); + lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult5 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0, quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0, quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0, quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult6 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.0, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1.0, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + { + matrixBuilder = storm::storage::SparseMatrixBuilder(22, 15, 28, true, true, 15); + matrixBuilder.newRowGroup(0); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.newRowGroup(1); + matrixBuilder.addNextValue(1, 0, 1); + matrixBuilder.addNextValue(2, 2, 1); + matrixBuilder.addNextValue(3, 4, 0.7); + matrixBuilder.addNextValue(3, 6, 0.3); + matrixBuilder.newRowGroup(4); + matrixBuilder.addNextValue(4, 0, 1); + + matrixBuilder.newRowGroup(5); + matrixBuilder.addNextValue(5, 4, 1); + matrixBuilder.addNextValue(6, 5, 0.8); + matrixBuilder.addNextValue(6, 9, 0.2); + matrixBuilder.newRowGroup(7); + matrixBuilder.addNextValue(7, 3, 1); + matrixBuilder.addNextValue(8, 5, 1); + matrixBuilder.newRowGroup(9); + matrixBuilder.addNextValue(9, 3, 1); + + matrixBuilder.newRowGroup(10); + matrixBuilder.addNextValue(10, 7, 1); + matrixBuilder.newRowGroup(11); + matrixBuilder.addNextValue(11, 6, 1); + matrixBuilder.addNextValue(12, 8, 1); + matrixBuilder.newRowGroup(13); + matrixBuilder.addNextValue(13, 6, 1); + + matrixBuilder.newRowGroup(14); + matrixBuilder.addNextValue(14, 10, 1); + matrixBuilder.newRowGroup(15); + matrixBuilder.addNextValue(15, 9, 1); + matrixBuilder.addNextValue(16, 11, 1); + matrixBuilder.newRowGroup(17); + matrixBuilder.addNextValue(17, 9, 1); + + matrixBuilder.newRowGroup(18); + matrixBuilder.addNextValue(18, 5, 0.4); + matrixBuilder.addNextValue(18, 8, 0.3); + matrixBuilder.addNextValue(18, 11, 0.3); + + matrixBuilder.newRowGroup(19); + matrixBuilder.addNextValue(19, 7, 0.7); + matrixBuilder.addNextValue(19, 12, 0.3); + + matrixBuilder.newRowGroup(20); + matrixBuilder.addNextValue(20, 12, 0.1); + matrixBuilder.addNextValue(20, 13, 0.9); + matrixBuilder.addNextValue(21, 12, 1); + + storm::storage::SparseMatrix 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::Mdp(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseMdpPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(storm::logic::OptimalityType::Maximize, labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(37./60., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(2./3., quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.5, quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1./3., quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(31./60., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(101./200., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(31./60., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision()); + + lraFormula = std::make_shared(storm::logic::OptimalityType::Minimize, labelFormula); + + result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1./3., quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.4, quantitativeResult2[3], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(1./3., quantitativeResult2[6], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[9], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.26, quantitativeResult2[12], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(467./1500., quantitativeResult2[13], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.26, quantitativeResult2[14], storm::settings::nativeEquationSolverSettings().getPrecision()); + } +} \ No newline at end of file From 5acaed6048e7ad832b7c107cd6e80f02b8e5885a Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 7 May 2015 13:07:00 +0200 Subject: [PATCH 17/20] Added flag to keep zeros when transposing. Former-commit-id: 811f6824cfebc45c1a5b9560bdc65a3201f3350b --- src/storage/SparseMatrix.cpp | 9 +++++---- src/storage/SparseMatrix.h | 7 ++++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 880d15d82..25beb4093 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -603,10 +603,11 @@ namespace storm { } template - SparseMatrix SparseMatrix::transpose(bool joinGroups) const { + SparseMatrix SparseMatrix::transpose(bool joinGroups, bool keepZeros) const { index_type rowCount = this->getColumnCount(); index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount(); - index_type entryCount = this->getEntryCount(); + //index_type entryCount = keepZeros ? this->getEntryCount() : this->getNonzeroEntryCount(); this wont work if someone modified the matrix after creation... + index_type entryCount = this->getEntryCount(); std::vector rowIndications(rowCount + 1); std::vector> columnsAndValues(entryCount); @@ -614,7 +615,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()) { + if (transition.getValue() != storm::utility::zero() || keepZeros) { ++rowIndications[transition.getColumn() + 1]; } } @@ -633,7 +634,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()) { + if (transition.getValue() != storm::utility::zero() || keepZeros) { columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); nextIndices[transition.getColumn()]++; } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 8eded966b..bf0193c7e 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -574,12 +574,13 @@ namespace storm { /*! * Transposes the matrix. - * - * @param joinGroups A flag indicating whether the row groups are supposed to be treated as single rows. + * + * @param joinGroups A flag indicating whether the row groups are supposed to be treated as single rows. + * @param keepZeros A flag indicating whether entries with value zero should be kept. * * @return A sparse matrix that represents the transpose of this matrix. */ - storm::storage::SparseMatrix transpose(bool joinGroups = false) const; + storm::storage::SparseMatrix transpose(bool joinGroups = false, bool keepZeros = false) const; /*! * Transforms the matrix into an equation system. That is, it transforms the matrix A into a matrix (1-A). From 8e688f71ff8dff97d60bac22dce91af99a029b6b Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 7 May 2015 13:07:40 +0200 Subject: [PATCH 18/20] Tests for DTMC LRA and some bugfixes. All tests pass. Former-commit-id: 589db6c2b3d1fbd6c34fef615b19bef686bee2f8 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 12 +- .../SparseDtmcPrctlModelCheckerTest.cpp | 145 ++++++++++++++++++ 2 files changed, 152 insertions(+), 5 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 9c70023a7..fdcf0bb3a 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -321,6 +321,8 @@ namespace storm { // Get some data members for convenience. typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); // First we check which states are in BSCCs. We use this later to speed up reachability analysis storm::storage::BitVector statesInBsccs(numOfStates); @@ -340,11 +342,9 @@ namespace storm { 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 bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true).transpose(); - - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); + storm::storage::SparseMatrix 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) { @@ -352,6 +352,8 @@ namespace storm { } } } + //now transpose, this internally removes all explicit zeros from the matrix that where introduced when subtracting the identity matrix + bsccEquationSystem = bsccEquationSystem.transpose(); std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); @@ -392,7 +394,7 @@ namespace storm { rewardRightSide.push_back(reward); } - storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs); + storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); rewardEquationSystemMatrix.convertToEquationSystem(); std::vector rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one); diff --git a/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp index c204d36dd..7f53f45b4 100644 --- a/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp @@ -127,3 +127,148 @@ TEST(SparseDtmcPrctlModelCheckerTest, SynchronousLeader) { EXPECT_NEAR(1.0448979589010925, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision()); } + +TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> dtmc; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(2, 2, 2); + matrixBuilder.addNextValue(0, 1, 1.); + matrixBuilder.addNextValue(1, 0, 1.); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + { + matrixBuilder = storm::storage::SparseMatrixBuilder(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 transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision()); + } + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(3, 3, 3); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 2, 1); + matrixBuilder.addNextValue(2, 0, 1); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 2); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + 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(SparseDtmcPrctlModelCheckerTest, LRA) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> mdp; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(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 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(transitionMatrix, ap, boost::none, boost::none, boost::none)); + + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + auto labelFormula = std::make_shared("a"); + auto lraFormula = std::make_shared(labelFormula); + + std::unique_ptr result = std::move(checker.check(*lraFormula)); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + 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()); + } +} \ No newline at end of file From cf5442fe4506383cb169cf89c479c30ed7e97af2 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 7 May 2015 16:23:29 +0200 Subject: [PATCH 19/20] Bugfix and test-fix: Only the "never leave MEC"-states have cost > 0 and transition costs are all 0 in the ssp. Former-commit-id: f6688a8956a17505554a3cc26e1e665165078547 --- .../csl/SparseMarkovAutomatonCslModelChecker.cpp | 9 +-------- src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp | 9 +-------- .../modelchecker/SparseMdpPrctlModelCheckerTest.cpp | 10 +++++----- 3 files changed, 7 insertions(+), 21 deletions(-) diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp index ebedb6c25..e3149901f 100644 --- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp @@ -354,14 +354,7 @@ namespace storm { // Now insert all (cumulative) probability values that target an MEC. for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) { if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) { - // If the target MEC is the same as the current one, instead of adding a transition, we need to add the weighted reward - // to the right-hand side vector of the SSP. - if (mecIndex == targetMecIndex) { - b.back() += auxiliaryStateToProbabilityMap[mecIndex] * lraValuesForEndComponents[mecIndex]; - } else { - // Otherwise, we add a transition to the auxiliary state that is associated with the target MEC. - sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]); - } + sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]); } } diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 00045105d..0a9007c22 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -440,14 +440,7 @@ namespace storm { // Now insert all (cumulative) probability values that target an MEC. for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) { if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) { - // If the target MEC is the same as the current one, instead of adding a transition, we need to add the weighted reward - // to the right-hand side vector of the SSP. - if (mecIndex == targetMecIndex) { - b.back() += auxiliaryStateToProbabilityMap[mecIndex] * lraValuesForEndComponents[mecIndex]; - } else { - // Otherwise, we add a transition to the auxiliary state that is associated with the target MEC. - sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]); - } + sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]); } } diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index d3210a395..59b2ad23f 100644 --- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp @@ -518,12 +518,12 @@ TEST(SparseMdpPrctlModelCheckerTest, LRA) { result = std::move(checker.check(*lraFormula)); storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); - EXPECT_NEAR(1./3., quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(0.4, quantitativeResult2[3], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult2[3], storm::settings::nativeEquationSolverSettings().getPrecision()); EXPECT_NEAR(1./3., quantitativeResult2[6], storm::settings::nativeEquationSolverSettings().getPrecision()); EXPECT_NEAR(0.0, quantitativeResult2[9], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(.26, quantitativeResult2[12], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(467./1500., quantitativeResult2[13], storm::settings::nativeEquationSolverSettings().getPrecision()); - EXPECT_NEAR(.26, quantitativeResult2[14], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult2[12], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(.79 / 3., quantitativeResult2[13], storm::settings::nativeEquationSolverSettings().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult2[14], storm::settings::nativeEquationSolverSettings().getPrecision()); } } \ No newline at end of file From 1cf0a73c4e009e103385492dd3a783c8e0c8cd31 Mon Sep 17 00:00:00 2001 From: David_Korzeniewski Date: Thu, 7 May 2015 18:18:26 +0200 Subject: [PATCH 20/20] Added methods to update nonzero entry count and update it when necessary And a fix for a compile error on gcc&clang. Former-commit-id: 2a095ca86493e5163d61b4a4262ff073ff1a7f38 --- .../prctl/SparseDtmcPrctlModelChecker.cpp | 2 +- src/storage/SparseMatrix.cpp | 55 +++++++++++++------ src/storage/SparseMatrix.h | 28 +++++++--- 3 files changed, 61 insertions(+), 24 deletions(-) diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index fdcf0bb3a..8f63110c8 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -451,7 +451,7 @@ namespace storm { } template - std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType = boost::optional()) { + std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 25beb4093..58a842cd3 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -239,20 +239,12 @@ namespace storm { template SparseMatrix::SparseMatrix(index_type columnCount, std::vector const& rowIndications, std::vector> const& columnsAndValues, std::vector const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) { - for (auto const& element : *this) { - if (element.getValue() != storm::utility::zero()) { - ++this->nonzeroEntryCount; - } - } + this->updateNonzeroEntryCount(); } template SparseMatrix::SparseMatrix(index_type columnCount, std::vector&& rowIndications, std::vector>&& columnsAndValues, std::vector&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) { - for (auto const& element : *this) { - if (element.getValue() != storm::utility::zero()) { - ++this->nonzeroEntryCount; - } - } + this->updateNonzeroEntryCount(); } template @@ -366,7 +358,23 @@ namespace storm { template typename SparseMatrix::index_type SparseMatrix::getNonzeroEntryCount() const { return nonzeroEntryCount; - } + } + + template + void SparseMatrix::updateNonzeroEntryCount() const { + //SparseMatrix* self = const_cast*>(this); + this->nonzeroEntryCount = 0; + for (auto const& element : *this) { + if (element.getValue() != storm::utility::zero()) { + ++this->nonzeroEntryCount; + } + } + } + + template + void SparseMatrix::updateNonzeroEntryCount(std::make_signed::type difference) { + this->nonzeroEntryCount += difference; + } template typename SparseMatrix::index_type SparseMatrix::getRowGroupCount() const { @@ -416,6 +424,7 @@ namespace storm { columnValuePtr->setValue(storm::utility::one()); ++columnValuePtr; for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) { + ++this->nonzeroEntryCount; columnValuePtr->setColumn(0); columnValuePtr->setValue(storm::utility::zero()); } @@ -606,8 +615,12 @@ namespace storm { SparseMatrix SparseMatrix::transpose(bool joinGroups, bool keepZeros) const { index_type rowCount = this->getColumnCount(); index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount(); - //index_type entryCount = keepZeros ? this->getEntryCount() : this->getNonzeroEntryCount(); this wont work if someone modified the matrix after creation... - index_type entryCount = this->getEntryCount(); + if (keepZeros) { + index_type entryCount = this->getEntryCount(); + } else { + this->updateNonzeroEntryCount(); + index_type entryCount = this->getNonzeroEntryCount(); + } std::vector rowIndications(rowCount + 1); std::vector> columnsAndValues(entryCount); @@ -660,13 +673,22 @@ namespace storm { template void SparseMatrix::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(); + // If there is a row without the diagonal element, an exception is thrown. + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); bool foundDiagonalElement = false; for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { - 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; } } @@ -696,6 +718,7 @@ namespace storm { for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { + --this->nonzeroEntryCount; entry.setValue(storm::utility::zero()); } } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index bf0193c7e..45e7b19ae 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -465,12 +465,26 @@ namespace storm { */ uint_fast64_t getRowGroupEntryCount(uint_fast64_t const group) const; - /*! - * Returns the number of nonzero entries in the matrix. - * - * @return The number of nonzero entries in the matrix. - */ - index_type getNonzeroEntryCount() const; + /*! + * Returns the cached number of nonzero entries in the matrix. + * + * @see updateNonzeroEntryCount() + * + * @return The number of nonzero entries in the matrix. + */ + index_type getNonzeroEntryCount() const; + + /*! + * Recompute the nonzero entry count + */ + void updateNonzeroEntryCount() const; + + /*! + * Change the nonzero entry count by the provided value. + * + * @param difference Difference between old and new nonzero entry count. + */ + void updateNonzeroEntryCount(std::make_signed::type difference); /*! * Returns the number of row groups in the matrix. @@ -827,7 +841,7 @@ namespace storm { index_type entryCount; // The number of nonzero entries in the matrix. - index_type nonzeroEntryCount; + mutable index_type nonzeroEntryCount; // The storage for the columns and values of all entries in the matrix. std::vector> columnsAndValues;