diff --git a/src/storm/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/storm/modelchecker/csl/HybridCtmcCslModelChecker.cpp index 88f380bd4..c4102efab 100644 --- a/src/storm/modelchecker/csl/HybridCtmcCslModelChecker.cpp +++ b/src/storm/modelchecker/csl/HybridCtmcCslModelChecker.cpp @@ -97,7 +97,7 @@ namespace storm { upperBound = storm::utility::infinity(); } - return storm::modelchecker::helper::HybridCtmcCslHelper::computeBoundedUntilProbabilities(env, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), checkTask.isQualitativeSet(), lowerBound, upperBound); + return storm::modelchecker::helper::HybridCtmcCslHelper::computeBoundedUntilProbabilities(env, this->getModel(), checkTask.isOnlyInitialStatesRelevantSet(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), checkTask.isQualitativeSet(), lowerBound, upperBound); } template @@ -105,7 +105,7 @@ namespace storm { storm::logic::InstantaneousRewardFormula const& rewardPathFormula = checkTask.getFormula(); STORM_LOG_THROW(!rewardPathFormula.isStepBounded(), storm::exceptions::NotImplementedException, "Currently step-bounded properties on CTMCs are not supported."); - return storm::modelchecker::helper::HybridCtmcCslHelper::computeInstantaneousRewards(env, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getRewardModel(""), rewardPathFormula.getBound()); + return storm::modelchecker::helper::HybridCtmcCslHelper::computeInstantaneousRewards(env, this->getModel(), checkTask.isOnlyInitialStatesRelevantSet(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getRewardModel(""), rewardPathFormula.getBound()); } template @@ -114,7 +114,7 @@ namespace storm { STORM_LOG_THROW(rewardPathFormula.getTimeBoundReference().isTimeBound(), storm::exceptions::NotImplementedException, "Currently step-bounded and reward-bounded properties on CTMCs are not supported."); auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask); - return storm::modelchecker::helper::HybridCtmcCslHelper::computeCumulativeRewards(env, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModel.get(), rewardPathFormula.getNonStrictBound()); + return storm::modelchecker::helper::HybridCtmcCslHelper::computeCumulativeRewards(env, this->getModel(), checkTask.isOnlyInitialStatesRelevantSet(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModel.get(), rewardPathFormula.getNonStrictBound()); } template @@ -123,13 +123,13 @@ namespace storm { std::unique_ptr subResultPointer = this->check(env, stateFormula); SymbolicQualitativeCheckResult const& subResult = subResultPointer->asSymbolicQualitativeCheckResult(); - return storm::modelchecker::helper::HybridCtmcCslHelper::computeLongRunAverageProbabilities(env, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), subResult.getTruthValuesVector()); + return storm::modelchecker::helper::HybridCtmcCslHelper::computeLongRunAverageProbabilities(env, this->getModel(), checkTask.isOnlyInitialStatesRelevantSet(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), subResult.getTruthValuesVector()); } template std::unique_ptr HybridCtmcCslModelChecker::computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) { auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask); - return storm::modelchecker::helper::HybridCtmcCslHelper::computeLongRunAverageRewards(env, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModel.get()); + return storm::modelchecker::helper::HybridCtmcCslHelper::computeLongRunAverageRewards(env, this->getModel(), checkTask.isOnlyInitialStatesRelevantSet(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModel.get()); } // Explicitly instantiate the model checker. diff --git a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp index 3555c7216..9d3972179 100644 --- a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp @@ -3,6 +3,8 @@ #include "storm/modelchecker/csl/helper/SparseCtmcCslHelper.h" #include "storm/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h" +#include "storm/environment/solver/TimeBoundedSolverEnvironment.h" + #include "storm/storage/dd/DdManager.h" #include "storm/storage/dd/Add.h" #include "storm/storage/dd/Bdd.h" @@ -45,13 +47,34 @@ namespace storm { } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound) { + std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound) { // If the time bounds are [0, inf], we rather call untimed reachability. if (storm::utility::isZero(lowerBound) && upperBound == storm::utility::infinity()) { return computeUntilProbabilities(env, model, rateMatrix, exitRateVector, phiStates, psiStates, qualitative); } + if (env.solver().isForceSoundness() && env.solver().timeBounded().getRelativeTerminationCriterion()) { + // Forward this query to the sparse engine + storm::utility::Stopwatch conversionWatch(true); + storm::dd::Odd odd = model.getReachableStates().createOdd(); + storm::storage::SparseMatrix explicitRateMatrix = rateMatrix.toMatrix(odd, odd); + std::vector explicitExitRateVector = exitRateVector.toVector(odd); + storm::solver::SolveGoal goal; + if (onlyInitialStatesRelevant) { + goal.setRelevantValues(model.getInitialStates().toVector(odd)); + } + conversionWatch.stop(); + STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); + + std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeBoundedUntilProbabilities(env, std::move(goal), explicitRateMatrix, explicitRateMatrix.transpose(true), phiStates.toVector(odd), psiStates.toVector(odd), explicitExitRateVector, qualitative, lowerBound, upperBound); + + return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); + } + + // Set the possible (absolute) error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()) / 8.0; + // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0 // or t' != inf. @@ -94,7 +117,7 @@ namespace storm { // Finally compute the transient probabilities. std::vector values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero()); - std::vector subresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values); + std::vector subresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, epsilon); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), (psiStates || !statesWithProbabilityGreater0) && model.getReachableStates(), @@ -141,7 +164,7 @@ namespace storm { STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); // Compute the transient probabilities. - result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result); + result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, epsilon); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), !relevantStates && model.getReachableStates(), model.getManager().template getAddZero(), relevantStates, odd, result)); } else { @@ -169,7 +192,7 @@ namespace storm { // Compute the transient probabilities. std::vector values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero()); - std::vector subResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values); + std::vector subResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, epsilon); // Transform the explicit result to a hybrid check result, so we can easily convert it to // a symbolic qualitative format. @@ -209,7 +232,7 @@ namespace storm { conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult); + newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, epsilon); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), !relevantStates && model.getReachableStates(), model.getManager().template getAddZero(), relevantStates, odd, newSubresult)); } else { @@ -236,7 +259,7 @@ namespace storm { conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult); + newSubresult = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, epsilon); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), !statesWithProbabilityGreater0 && model.getReachableStates(), model.getManager().template getAddZero(), statesWithProbabilityGreater0, odd, newSubresult)); } @@ -248,12 +271,12 @@ namespace storm { } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound) { + std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound) { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded until probabilities is unsupported for this value type."); } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { + std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { // Only compute the result if the model has a state-based reward model. STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -266,10 +289,12 @@ namespace storm { conversionWatch.stop(); // Initialize result to state rewards of the model. - std::vector result = rewardModel.getStateRewardVector().toVector(odd); + auto rewardsAdd = rewardModel.getStateRewardVector(); + std::vector result = rewardsAdd.toVector(odd); + ValueType maxValue = std::max(rewardsAdd.getMax(), -rewardsAdd.getMin()); - // If the time-bound is not zero, we need to perform a transient analysis. - if (timeBound > 0) { + // If the rewards are not zero and the time-bound is not zero, we need to perform a transient analysis. + if (!storm::utility::isZero(maxValue) && timeBound > 0) { ValueType uniformizationRate = 1.02 * exitRateVector.getMax(); STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); @@ -280,19 +305,39 @@ namespace storm { conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result); + // Set the possible error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + if (env.solver().timeBounded().getRelativeTerminationCriterion()) { + // Be more precise, if the maximum value is very small (This still gives no sound guarantee!) + epsilon *= std::min(storm::utility::one(), maxValue); + } else { + // Be more precise, if the maximal possible value is very large + epsilon /= std::max(storm::utility::one(), maxValue); + } + + storm::storage::BitVector relevantValues; + if (onlyInitialStatesRelevant) { + relevantValues = model.getInitialStates().toVector(odd); + } else { + relevantValues = storm::storage::BitVector(result.size(), true); + } + + // Loop until the desired precision is reached. + do { + result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, epsilon); + } while (storm::modelchecker::helper::SparseCtmcCslHelper::checkAndUpdateTransientProbabilityEpsilon(env, epsilon, result, relevantValues)); } return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), odd, result)); } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { + std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing instantaneous rewards is unsupported for this value type."); } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { + std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { // Only compute the result if the model has a state-based reward model. STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -327,18 +372,47 @@ namespace storm { conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); + ValueType maxReward = std::max(totalRewardVector.getMax(), -totalRewardVector.getMin()); + + // If all rewards are zero, the result is the constant zero vector. + if (storm::utility::isZero(maxReward)) { + return std::unique_ptr(new SymbolicQuantitativeCheckResult(model.getReachableStates(), model.getManager().template getAddZero())); + } + + // Set the possible (absolute) error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + if (env.solver().timeBounded().getRelativeTerminationCriterion()) { + // Be more precise, if the value is very small (this still gives no sound guarantee) + epsilon *= std::min(storm::utility::one(), maxReward); + } else { + // Be more precise, if the maximal possible value is very large + epsilon /= std::max(storm::utility::one(), maxReward * timeBound); + } + + storm::storage::BitVector relevantValues; + if (onlyInitialStatesRelevant) { + relevantValues = model.getInitialStates().toVector(odd); + } else { + relevantValues = storm::storage::BitVector(explicitTotalRewardVector.size(), true); + } + + std::vector result; // Finally, compute the transient probabilities. - std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector); + // Loop until the desired precision is reached. + do { + result = storm::modelchecker::helper::SparseCtmcCslHelper::computeTransientProbabilities(env, explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, epsilon); + } while (storm::modelchecker::helper::SparseCtmcCslHelper::checkAndUpdateTransientProbabilityEpsilon(env, epsilon, result, relevantValues)); + return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); } template::SupportsExponential, int>::type> - std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { + std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound) { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing cumulative rewards is unsupported for this value type."); } template - std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates) { + std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates) { storm::utility::Stopwatch conversionWatch(true); @@ -347,16 +421,20 @@ namespace storm { storm::storage::SparseMatrix explicitRateMatrix = rateMatrix.toMatrix(odd, odd); std::vector explicitExitRateVector = exitRateVector.toVector(odd); + storm::solver::SolveGoal goal; + if (onlyInitialStatesRelevant) { + goal.setRelevantValues(model.getInitialStates().toVector(odd)); + } conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal(), explicitRateMatrix, psiStates.toVector(odd), &explicitExitRateVector); + std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, std::move(goal), explicitRateMatrix, psiStates.toVector(odd), &explicitExitRateVector); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); } template - std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel) { + std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel) { STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); storm::dd::Add probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector); @@ -368,10 +446,14 @@ namespace storm { storm::storage::SparseMatrix explicitRateMatrix = rateMatrix.toMatrix(odd, odd); std::vector explicitExitRateVector = exitRateVector.toVector(odd); + storm::solver::SolveGoal goal; + if (onlyInitialStatesRelevant) { + goal.setRelevantValues(model.getInitialStates().toVector(odd)); + } conversionWatch.stop(); STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms."); - std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal(), explicitRateMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector); + std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, std::move(goal), explicitRateMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector); return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), std::move(odd), std::move(result))); } @@ -402,50 +484,50 @@ namespace storm { // Explicit instantiations. // Cudd, double. - template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); - template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); - template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template std::unique_ptr HybridCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative); template std::unique_ptr HybridCtmcCslHelper::computeReachabilityRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); template std::unique_ptr HybridCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& nextStates); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); template storm::dd::Add HybridCtmcCslHelper::computeProbabilityMatrix(storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector); template storm::dd::Add HybridCtmcCslHelper::computeUniformizedMatrix(storm::models::symbolic::Ctmc const& model, storm::dd::Add const& transitionMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& maybeStates, double uniformizationRate); // Sylvan, double. - template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); - template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); - template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template std::unique_ptr HybridCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative); template std::unique_ptr HybridCtmcCslHelper::computeReachabilityRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); template std::unique_ptr HybridCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& nextStates); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); template storm::dd::Add HybridCtmcCslHelper::computeProbabilityMatrix(storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector); template storm::dd::Add HybridCtmcCslHelper::computeUniformizedMatrix(storm::models::symbolic::Ctmc const& model, storm::dd::Add const& transitionMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& maybeStates, double uniformizationRate); // Sylvan, rational number. - template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); - template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); - template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template std::unique_ptr HybridCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative); template std::unique_ptr HybridCtmcCslHelper::computeReachabilityRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); template std::unique_ptr HybridCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& nextStates); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); template storm::dd::Add HybridCtmcCslHelper::computeProbabilityMatrix(storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector); template storm::dd::Add HybridCtmcCslHelper::computeUniformizedMatrix(storm::models::symbolic::Ctmc const& model, storm::dd::Add const& transitionMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& maybeStates, storm::RationalNumber uniformizationRate); // Sylvan, rational function. - template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); - template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); - template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + template std::unique_ptr HybridCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + template std::unique_ptr HybridCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template std::unique_ptr HybridCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative); template std::unique_ptr HybridCtmcCslHelper::computeReachabilityRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); template std::unique_ptr HybridCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& nextStates); - template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); + template std::unique_ptr HybridCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); template storm::dd::Add HybridCtmcCslHelper::computeProbabilityMatrix(storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector); template storm::dd::Add HybridCtmcCslHelper::computeUniformizedMatrix(storm::models::symbolic::Ctmc const& model, storm::dd::Add const& transitionMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& maybeStates, storm::RationalFunction uniformizationRate); diff --git a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.h b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.h index 8814ca9a7..73c5e3c46 100644 --- a/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.h +++ b/src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.h @@ -21,22 +21,22 @@ namespace storm { class HybridCtmcCslHelper { public: template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + static std::unique_ptr computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); + static std::unique_ptr computeBoundedUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, double lowerBound, double upperBound); template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + static std::unique_ptr computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + static std::unique_ptr computeInstantaneousRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + static std::unique_ptr computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template::SupportsExponential, int>::type = 0> - static std::unique_ptr computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); + static std::unique_ptr computeCumulativeRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, double timeBound); template static std::unique_ptr computeUntilProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative); @@ -45,13 +45,13 @@ namespace storm { static std::unique_ptr computeReachabilityRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel, storm::dd::Bdd const& targetStates, bool qualitative); template - static std::unique_ptr computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); + static std::unique_ptr computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& psiStates); template static std::unique_ptr computeNextProbabilities(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, storm::dd::Bdd const& nextStates); template - static std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); + static std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::models::symbolic::Ctmc const& model, bool onlyInitialStatesRelevant, storm::dd::Add const& rateMatrix, storm::dd::Add const& exitRateVector, typename storm::models::symbolic::Model::RewardModelType const& rewardModel); /*! * Converts the given rate-matrix into a time-abstract probability matrix. diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index ae8dc0bfc..20993a059 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -34,6 +34,39 @@ namespace storm { namespace modelchecker { namespace helper { + + template + bool SparseCtmcCslHelper::checkAndUpdateTransientProbabilityEpsilon(storm::Environment const& env, ValueType& epsilon, std::vector const& resultVector, storm::storage::BitVector const& relevantPositions) { + // Check if the check is necessary for the provided settings + if (!env.solver().isForceSoundness() || !env.solver().timeBounded().getRelativeTerminationCriterion()) { + // No need to update epsilon + return false; + } + + ValueType precision = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + // If we need to compute values with relative precision, it might be necessary to increase the precision requirements (epsilon) + ValueType newEpsilon = epsilon; + // Only consider positions that are relevant for the solve goal (e.g. initial states of the model) and are supposed to have a non-zero value + for (auto const& state : relevantPositions) { + if (storm::utility::isZero(resultVector[state])) { + newEpsilon = std::min(epsilon * storm::utility::convertNumber(0.1), newEpsilon); + } else { + ValueType relativeError = epsilon / resultVector[state]; // epsilon is an upper bound for the absolute error we made + if (relativeError > precision) { + newEpsilon = std::min(resultVector[state] * precision, newEpsilon); + } + } + } + if (newEpsilon < epsilon) { + STORM_LOG_INFO("Re-computing transient probabilities with new truncation error " << newEpsilon << " to guarantee sound results with relative precision."); + epsilon = newEpsilon; + return true; + } else { + return false; + } + } + + template ::SupportsExponential, int>::type> std::vector SparseCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound) { @@ -52,6 +85,9 @@ namespace storm { // Create the result vector. std::vector result; + // Set the possible (absolute) error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()) / 8.0; + // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the // further computations. storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(backwardTransitions, phiStates, psiStates); @@ -59,89 +95,38 @@ namespace storm { storm::storage::BitVector statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 & ~psiStates; STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNumberOfSetBits() << " 'maybe' states."); - if (!statesWithProbabilityGreater0.empty()) { - if (storm::utility::isZero(upperBound)) { - // In this case, the interval is of the form [0, 0]. - result = std::vector(numberOfStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); - } else { - if (storm::utility::isZero(lowerBound)) { - // In this case, the interval is of the form [0, t]. - // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier. - + // the positions within the result for which the precision needs to be checked + storm::storage::BitVector relevantValues; + if (goal.hasRelevantValues()) { + relevantValues = std::move(goal.relevantValues()); + relevantValues &= statesWithProbabilityGreater0; + } else { + relevantValues = statesWithProbabilityGreater0; + } + + do { // Iterate until the desired precision is reached (only relevant for relative precision criterion) + if (!statesWithProbabilityGreater0.empty()) { + if (storm::utility::isZero(upperBound)) { + // In this case, the interval is of the form [0, 0]. result = std::vector(numberOfStates, storm::utility::zero()); storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); - if (!statesWithProbabilityGreater0NonPsi.empty()) { - // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. - ValueType uniformizationRate = 0; - for (auto const& state : statesWithProbabilityGreater0NonPsi) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - // Compute the uniformized matrix. - storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); - - // Compute the vector that is to be added as a compensation for removing the absorbing states. - std::vector b = rateMatrix.getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); - for (auto& element : b) { - element /= uniformizationRate; - } - - // Finally compute the transient probabilities. - std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); - std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, &b, upperBound, uniformizationRate, values); - storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0NonPsi, subresult); - } - } else if (upperBound == storm::utility::infinity()) { - // In this case, the interval is of the form [t, inf] with t != 0. - - // Start by computing the (unbounded) reachability probabilities of reaching psi states while - // staying in phi states. - result = computeUntilProbabilities(env, storm::solver::SolveGoal(), rateMatrix, backwardTransitions, exitRates, phiStates, psiStates, qualitative); - - // Determine the set of states that must be considered further. - storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; - std::vector subResult(relevantStates.getNumberOfSetBits()); - storm::utility::vector::selectVectorValues(subResult, relevantStates, result); - - ValueType uniformizationRate = 0; - for (auto const& state : relevantStates) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - // Compute the uniformized matrix. - storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, relevantStates, uniformizationRate, exitRates); - - // Compute the transient probabilities. - subResult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, subResult); - - // Fill in the correct values. - storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, relevantStates, subResult); } else { - // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf. - - if (lowerBound != upperBound) { - // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'. + if (storm::utility::isZero(lowerBound)) { + // In this case, the interval is of the form [0, t]. + // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier. - - storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; - std::vector newSubresult(relevantStates.getNumberOfSetBits(), storm::utility::zero()); - storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one()); + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); if (!statesWithProbabilityGreater0NonPsi.empty()) { // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. - ValueType uniformizationRate = storm::utility::zero(); + ValueType uniformizationRate = 0; for (auto const& state : statesWithProbabilityGreater0NonPsi) { uniformizationRate = std::max(uniformizationRate, exitRates[state]); } uniformizationRate *= 1.02; STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - // Compute the (first) uniformized matrix. + // Compute the uniformized matrix. storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); // Compute the vector that is to be added as a compensation for removing the absorbing states. @@ -150,59 +135,121 @@ namespace storm { element /= uniformizationRate; } - // Start by computing the transient probabilities of reaching a psi state in time t' - t. + // Finally compute the transient probabilities. std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); - std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values); - storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); + std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, &b, upperBound, uniformizationRate, values, epsilon); + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0NonPsi, subresult); } + } else if (upperBound == storm::utility::infinity()) { + // In this case, the interval is of the form [t, inf] with t != 0. + + // Start by computing the (unbounded) reachability probabilities of reaching psi states while + // staying in phi states. + result = computeUntilProbabilities(env, storm::solver::SolveGoal(), rateMatrix, backwardTransitions, exitRates, phiStates, psiStates, qualitative); - // Then compute the transient probabilities of being in such a state after t time units. For this, - // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. - ValueType uniformizationRate = storm::utility::zero(); + // Determine the set of states that must be considered further. + storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; + std::vector subResult(relevantStates.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subResult, relevantStates, result); + + ValueType uniformizationRate = 0; for (auto const& state : relevantStates) { uniformizationRate = std::max(uniformizationRate, exitRates[state]); } uniformizationRate *= 1.02; STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - // Finally, we compute the second set of transient probabilities. + // Compute the uniformized matrix. storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, relevantStates, uniformizationRate, exitRates); - newSubresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult); + + // Compute the transient probabilities. + subResult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, subResult, epsilon); // Fill in the correct values. - result = std::vector(numberOfStates, storm::utility::zero()); storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, relevantStates, newSubresult); + storm::utility::vector::setVectorValues(result, relevantStates, subResult); } else { - // In this case, the interval is of the form [t, t] with t != 0, t != inf. - - std::vector newSubresult = std::vector(statesWithProbabilityGreater0.getNumberOfSetBits()); - storm::utility::vector::setVectorValues(newSubresult, psiStates % statesWithProbabilityGreater0, storm::utility::one()); + // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf. - // Then compute the transient probabilities of being in such a state after t time units. For this, - // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. - ValueType uniformizationRate = storm::utility::zero(); - for (auto const& state : statesWithProbabilityGreater0) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); + if (lowerBound != upperBound) { + // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'. + + + storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; + std::vector newSubresult(relevantStates.getNumberOfSetBits(), storm::utility::zero()); + storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one()); + if (!statesWithProbabilityGreater0NonPsi.empty()) { + // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. + ValueType uniformizationRate = storm::utility::zero(); + for (auto const& state : statesWithProbabilityGreater0NonPsi) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Compute the (first) uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); + + // Compute the vector that is to be added as a compensation for removing the absorbing states. + std::vector b = rateMatrix.getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); + for (auto& element : b) { + element /= uniformizationRate; + } + + // Start by computing the transient probabilities of reaching a psi state in time t' - t. + std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); + // divide the possible error by two since we will make this error two times. + std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values, epsilon / storm::utility::convertNumber(2.0)); + storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); + } + + // Then compute the transient probabilities of being in such a state after t time units. For this, + // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. + ValueType uniformizationRate = storm::utility::zero(); + for (auto const& state : relevantStates) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Finally, we compute the second set of transient probabilities. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, relevantStates, uniformizationRate, exitRates); + newSubresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, epsilon / storm::utility::convertNumber(2.0)); + + // Fill in the correct values. + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, relevantStates, newSubresult); + } else { + // In this case, the interval is of the form [t, t] with t != 0, t != inf. + + std::vector newSubresult = std::vector(statesWithProbabilityGreater0.getNumberOfSetBits()); + storm::utility::vector::setVectorValues(newSubresult, psiStates % statesWithProbabilityGreater0, storm::utility::one()); + + // Then compute the transient probabilities of being in such a state after t time units. For this, + // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. + ValueType uniformizationRate = storm::utility::zero(); + for (auto const& state : statesWithProbabilityGreater0) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Finally, we compute the second set of transient probabilities. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0, uniformizationRate, exitRates); + newSubresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, epsilon); + + // Fill in the correct values. + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, ~statesWithProbabilityGreater0, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, newSubresult); } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - // Finally, we compute the second set of transient probabilities. - storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0, uniformizationRate, exitRates); - newSubresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult); - - // Fill in the correct values. - result = std::vector(numberOfStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, ~statesWithProbabilityGreater0, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, newSubresult); } } + } else { + result = std::vector(numberOfStates, storm::utility::zero()); } - } else { - result = std::vector(numberOfStates, storm::utility::zero()); - } - + } while (checkAndUpdateTransientProbabilityEpsilon(env, epsilon, result, relevantValues)); return result; } @@ -236,19 +283,48 @@ namespace storm { // Initialize result to state rewards of the this->getModel(). std::vector result(rewardModel.getStateRewardVector()); - // If the time-bound is not zero, we need to perform a transient analysis. - if (timeBound > 0) { - ValueType uniformizationRate = 0; - for (auto const& rate : exitRateVector) { - uniformizationRate = std::max(uniformizationRate, rate); - } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, storm::storage::BitVector(numberOfStates, true), uniformizationRate, exitRateVector); - result = computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, result); + // If the time-bound is zero, just return the current reward vector + if (storm::utility::isZero(timeBound)) { + return result; } + ValueType maxValue = storm::utility::vector::maximumElementAbs(result); + // If all entries are zero, we return the zero-vector + if (storm::utility::isZero(maxValue)) { + return result; + } + + ValueType uniformizationRate = 0; + for (auto const& rate : exitRateVector) { + uniformizationRate = std::max(uniformizationRate, rate); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, storm::storage::BitVector(numberOfStates, true), uniformizationRate, exitRateVector); + + // Set the possible error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + if (env.solver().timeBounded().getRelativeTerminationCriterion()) { + // Be more precise, if the maximum value is very small (precision can/has to be refined later) + epsilon *= std::min(storm::utility::one(), maxValue); + } else { + // Be more precise, if the maximal possible value is very large + epsilon /= std::max(storm::utility::one(), maxValue); + } + + storm::storage::BitVector relevantValues; + if (goal.hasRelevantValues()) { + relevantValues = std::move(goal.relevantValues()); + } else { + relevantValues = storm::storage::BitVector(result.size(), true); + } + + // Loop until the desired precision is reached. + do { + result = computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, result, epsilon); + } while (checkAndUpdateTransientProbabilityEpsilon(env, epsilon, result, relevantValues)); + return result; } @@ -285,9 +361,37 @@ namespace storm { // Compute the total state reward vector. std::vector totalRewardVector = rewardModel.getTotalRewardVector(rateMatrix, exitRateVector); + ValueType maxReward = storm::utility::vector::maximumElementAbs(totalRewardVector); + + // If all rewards are zero, the result is the constant zero vector. + if (storm::utility::isZero(maxReward)) { + return std::vector(numberOfStates, storm::utility::zero()); + } + + // Set the possible (absolute) error allowed for truncation (epsilon for fox-glynn) + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + if (env.solver().timeBounded().getRelativeTerminationCriterion()) { + // Be more precise, if the value is very small (precision can/has to be refined later) + epsilon *= std::min(storm::utility::one(), maxReward); + } else { + // Be more precise, if the maximal possible value is very large + epsilon /= std::max(storm::utility::one(), maxReward * timeBound); + } + + storm::storage::BitVector relevantValues; + if (goal.hasRelevantValues()) { + relevantValues = std::move(goal.relevantValues()); + } else { + relevantValues = storm::storage::BitVector(totalRewardVector.size(), true); + } - // Finally, compute the transient probabilities. - return computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, totalRewardVector); + // Loop until the desired precision is reached. + std::vector result; + do { + result = computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, totalRewardVector, epsilon); + } while (checkAndUpdateTransientProbabilityEpsilon(env, epsilon, result, relevantValues)); + + return result; } template ::SupportsExponential, int>::type> @@ -972,7 +1076,9 @@ namespace storm { element /= uniformizationRate; std::cout << element << std::endl; }*/ - + + ValueType epsilon = storm::utility::convertNumber(env.solver().timeBounded().getPrecision()) / 8.0; + STORM_LOG_WARN_COND(!env.solver().timeBounded().getRelativeTerminationCriterion(), "Computation of transient probabilities with relative precision not supported. Using absolute precision instead."); std::vector values(relevantStates.getNumberOfSetBits(), storm::utility::zero()); // Set initial states size_t i = 0; @@ -984,7 +1090,7 @@ namespace storm { ++i; } // Finally compute the transient probabilities. - std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, values); + std::vector subresult = computeTransientProbabilities(env, uniformizedMatrix, nullptr, timeBound, uniformizationRate, values, epsilon); storm::utility::vector::setVectorValues(result, relevantStates, subresult); } @@ -1025,8 +1131,8 @@ namespace storm { } template::SupportsExponential, int>::type> - std::vector SparseCtmcCslHelper::computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values) { - + std::vector SparseCtmcCslHelper::computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, ValueType epsilon) { + STORM_LOG_WARN_COND(epsilon > storm::utility::convertNumber(1e-20), "Very low truncation error " << epsilon << " requested. Numerical inaccuracies are possible."); ValueType lambda = timeBound * uniformizationRate; // If no time can pass, the current values are the result. @@ -1035,9 +1141,7 @@ namespace storm { } // Use Fox-Glynn to get the truncation points and the weights. -// std::tuple> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule().getPrecision() / 8.0); - - storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::utility::convertNumber(env.solver().timeBounded().getPrecision()) / 8.0); + storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, epsilon); STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << foxGlynnResult.left << ", right=" << foxGlynnResult.right); // Scale the weights so they add up to one. @@ -1157,7 +1261,7 @@ namespace storm { template storm::storage::SparseMatrix SparseCtmcCslHelper::computeUniformizedMatrix(storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& maybeStates, double uniformizationRate, std::vector const& exitRates); - template std::vector SparseCtmcCslHelper::computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, double timeBound, double uniformizationRate, std::vector values); + template std::vector SparseCtmcCslHelper::computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, double timeBound, double uniformizationRate, std::vector values, double epsilon); #ifdef STORM_HAVE_CARL template std::vector SparseCtmcCslHelper::computeBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound); diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h index 95840ae19..697eca743 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h @@ -93,13 +93,13 @@ namespace storm { * @param timeBound The time bound to use. * @param uniformizationRate The used uniformization rate. * @param values A vector mapping each state to an initial probability. - * @param linearEquationSolverFactory The factory to use when instantiating new linear equation solvers. - * @param useMixedPoissonProbabilities If set to true, instead of taking the poisson probabilities, mixed + * @param epsilon The precision used for computing the truncation points + * @tparam useMixedPoissonProbabilities If set to true, instead of taking the poisson probabilities, mixed * poisson probabilities are used. * @return The vector of transient probabilities. */ template::SupportsExponential, int>::type = 0> - static std::vector computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values); + static std::vector computeTransientProbabilities(Environment const& env, storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, ValueType epsilon); /*! * Converts the given rate-matrix into a time-abstract probability matrix. @@ -121,6 +121,16 @@ namespace storm { template static storm::storage::SparseMatrix computeGeneratorMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates); + /*! + * Checks whether the given result vector is sufficiently precise, according to the provided epsilon and the specified settings + * @param epsilon The truncation error. If the result needs to be more precise, this value will be decreased + * @param resultVector the currently obtained result + * @param relevantPositions Only these positions of the resultVector will be considered. + * @return + */ + template + static bool checkAndUpdateTransientProbabilityEpsilon(storm::Environment const& env, ValueType& epsilon, std::vector const& resultVector, storm::storage::BitVector const& relevantPositions); + private: template static std::vector computeLongRunAverages(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& rateMatrix, std::function const& valueGetter, std::vector const* exitRateVector); diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 0deed45ba..e5fda444e 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -928,6 +928,16 @@ namespace storm { return true; } + + template + T maximumElementAbs(std::vector const& vector) { + T res = storm::utility::zero(); + for (auto const& element : vector) { + res = std::max(res, storm::utility::abs(element)); + } + return res; + } + template T maximumElementDiff(std::vector const& vectorLeft, std::vector const& vectorRight) { T maxDiff = storm::utility::zero();