From 97936cbd8eb214f62059aab0f7f9277f85efd083 Mon Sep 17 00:00:00 2001 From: masawei Date: Thu, 26 Feb 2015 13:25:18 +0100 Subject: [PATCH 1/9] Found a fix for a bug causing the functional tests to segfault at DeterministicModelBisimulationDecomposition.Die. - By setting the blocks to be not sorted and unique a different constructor is used by the boost container. This prevents the segfault. |- I can't say exactly why this works nor do I know if the blocks are actually sorted and unique in the sense meant by the underlying container implementation. Former-commit-id: a1bfbab75a24c90cd0a54ac1e8bec5cabb9fe978 --- src/storage/DeterministicModelBisimulationDecomposition.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/storage/DeterministicModelBisimulationDecomposition.cpp b/src/storage/DeterministicModelBisimulationDecomposition.cpp index 6e974f58c..b20311ba9 100644 --- a/src/storage/DeterministicModelBisimulationDecomposition.cpp +++ b/src/storage/DeterministicModelBisimulationDecomposition.cpp @@ -875,7 +875,7 @@ namespace storm { // Convert the state-value-pairs to states only. std::function const&)> projection = [] (std::pair const& a) -> storm::storage::sparse::state_type { return a.first; }; - this->blocks[block.getId()] = block_type(boost::make_transform_iterator(partition.getBegin(block), projection), boost::make_transform_iterator(partition.getEnd(block), projection), true); + this->blocks[block.getId()] = block_type(boost::make_transform_iterator(partition.getBegin(block), projection), boost::make_transform_iterator(partition.getEnd(block), projection), false); } // If we are required to build the quotient model, do so now. @@ -1493,4 +1493,4 @@ namespace storm { template class DeterministicModelBisimulationDecomposition; #endif } -} \ No newline at end of file +} From 7e14dc031b89364a232231be6c88f3c33c62d7a2 Mon Sep 17 00:00:00 2001 From: dehnert Date: Thu, 26 Feb 2015 13:35:10 +0100 Subject: [PATCH 2/9] Reverted the last commit. The flag is there for performance reasons and there is no reason why it shouldn't work that way. Former-commit-id: e551eb461f3952a7f94e8fb7f27221526f5c35d3 --- src/storage/DeterministicModelBisimulationDecomposition.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/storage/DeterministicModelBisimulationDecomposition.cpp b/src/storage/DeterministicModelBisimulationDecomposition.cpp index b20311ba9..58d7c618e 100644 --- a/src/storage/DeterministicModelBisimulationDecomposition.cpp +++ b/src/storage/DeterministicModelBisimulationDecomposition.cpp @@ -875,7 +875,7 @@ namespace storm { // Convert the state-value-pairs to states only. std::function const&)> projection = [] (std::pair const& a) -> storm::storage::sparse::state_type { return a.first; }; - this->blocks[block.getId()] = block_type(boost::make_transform_iterator(partition.getBegin(block), projection), boost::make_transform_iterator(partition.getEnd(block), projection), false); + this->blocks[block.getId()] = block_type(boost::make_transform_iterator(partition.getBegin(block), projection), boost::make_transform_iterator(partition.getEnd(block), projection), true); } // If we are required to build the quotient model, do so now. From 546e047b8deb00769d921dcb3ea671e033b21c38 Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 16 Mar 2015 17:40:11 +0100 Subject: [PATCH 3/9] Fixed a bug that prevented correct comparison with bounds in formulas. Former-commit-id: ae6c28dcbe8057f8b087bedf06151fcebbada3ed --- .../results/ExplicitQuantitativeCheckResult.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp b/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp index 31b7e5296..3cb22b151 100644 --- a/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp +++ b/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp @@ -127,28 +127,28 @@ namespace storm { switch (comparisonType) { case logic::Less: for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) { - if (result[index] < bound) { + if (valuesAsVector[index] < bound) { result.set(index); } } break; case logic::LessEqual: for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) { - if (result[index] <= bound) { + if (valuesAsVector[index] <= bound) { result.set(index); } } break; case logic::Greater: for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) { - if (result[index] > bound) { + if (valuesAsVector[index] > bound) { result.set(index); } } break; case logic::GreaterEqual: for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) { - if (result[index] >= bound) { + if (valuesAsVector[index] >= bound) { result.set(index); } } From a44a3554c8053f288350b99eed73959ced3bc856 Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 16 Mar 2015 22:47:09 +0100 Subject: [PATCH 4/9] Fixed minimal command counterexample generation. Former-commit-id: 6e7e6208daacd31321d29ec3eec25cee56bc62e6 --- .../MILPMinimalLabelSetGenerator.h | 9 +++--- .../SMTMinimalCommandSetGenerator.h | 5 ++- src/modelchecker/AbstractModelChecker.cpp | 32 ++++++++++++++++--- .../prctl/SparseMdpPrctlModelChecker.cpp | 12 +++---- .../prctl/SparseMdpPrctlModelChecker.h | 6 +++- src/utility/cli.h | 5 +++ 6 files changed, 51 insertions(+), 18 deletions(-) diff --git a/src/counterexamples/MILPMinimalLabelSetGenerator.h b/src/counterexamples/MILPMinimalLabelSetGenerator.h index 4e09ae838..4c9b5b566 100644 --- a/src/counterexamples/MILPMinimalLabelSetGenerator.h +++ b/src/counterexamples/MILPMinimalLabelSetGenerator.h @@ -629,7 +629,7 @@ namespace storm { uint_fast64_t numberOfConstraintsCreated = 0; for (auto label : choiceInformation.knownLabels) { - storm::expressions::Expression constraint = variableInformation.labelToVariableMap.at(label) == solver.getConstant(0); + storm::expressions::Expression constraint = variableInformation.labelToVariableMap.at(label) == solver.getConstant(1); solver.addConstraint("KnownLabels" + std::to_string(numberOfConstraintsCreated), constraint); ++numberOfConstraintsCreated; } @@ -929,10 +929,9 @@ namespace storm { double maximalReachabilityProbability = 0; if (checkThresholdFeasible) { storm::modelchecker::SparseMdpPrctlModelChecker modelchecker(labeledMdp); - std::unique_ptr result = modelchecker.check(pathFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult const& quantitativeResult = result->asExplicitQuantitativeCheckResult(); + std::vector result = modelchecker.computeUntilProbabilitiesHelper(false, phiStates, psiStates, false); for (auto state : labeledMdp.getInitialStates()) { - maximalReachabilityProbability = std::max(maximalReachabilityProbability, quantitativeResult[state]); + maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]); } STORM_LOG_THROW((strictBound && maximalReachabilityProbability >= probabilityThreshold) || (!strictBound && maximalReachabilityProbability > probabilityThreshold), storm::exceptions::InvalidArgumentException, "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << "."); std::cout << std::endl << "Maximal reachability in model is " << maximalReachabilityProbability << "." << std::endl << std::endl; @@ -953,6 +952,8 @@ namespace storm { // (4.2) Construct constraint system. buildConstraintSystem(*solver, labeledMdp, psiStates, stateInformation, choiceInformation, variableInformation, probabilityThreshold, strictBound, includeSchedulerCuts); + solver->writeModelToFile("model.lp"); + // (4.3) Optimize the model. solver->optimize(); diff --git a/src/counterexamples/SMTMinimalCommandSetGenerator.h b/src/counterexamples/SMTMinimalCommandSetGenerator.h index bde97a841..e79d76dda 100644 --- a/src/counterexamples/SMTMinimalCommandSetGenerator.h +++ b/src/counterexamples/SMTMinimalCommandSetGenerator.h @@ -1627,10 +1627,9 @@ namespace storm { double maximalReachabilityProbability = 0; if (checkThresholdFeasible) { storm::modelchecker::SparseMdpPrctlModelChecker modelchecker(labeledMdp); - std::unique_ptr result = modelchecker.check(pathFormula); - storm::modelchecker::ExplicitQuantitativeCheckResult const& quantitativeResult = result->asExplicitQuantitativeCheckResult(); + std::vector result = modelchecker.computeUntilProbabilitiesHelper(false, phiStates, psiStates, false); for (auto state : labeledMdp.getInitialStates()) { - maximalReachabilityProbability = std::max(maximalReachabilityProbability, quantitativeResult[state]); + maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]); } STORM_LOG_THROW((strictBound && maximalReachabilityProbability >= probabilityThreshold) || (!strictBound && maximalReachabilityProbability > probabilityThreshold), storm::exceptions::InvalidArgumentException, "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << "."); std::cout << std::endl << "Maximal reachability in model is " << maximalReachabilityProbability << "." << std::endl << std::endl; diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp index eb9d539a9..6b4ab55ab 100644 --- a/src/modelchecker/AbstractModelChecker.cpp +++ b/src/modelchecker/AbstractModelChecker.cpp @@ -85,7 +85,7 @@ namespace storm { std::unique_ptr AbstractModelChecker::computeLongRunAverage(storm::logic::StateFormula const& eventuallyFormula, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker does not support the computation of long-run averages."); } - + std::unique_ptr AbstractModelChecker::computeExpectedTimes(storm::logic::EventuallyFormula const& eventuallyFormula, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker does not support the computation of expected times."); } @@ -160,6 +160,12 @@ namespace storm { std::unique_ptr result; if (stateFormula.hasOptimalityType()) { result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, stateFormula.getOptimalityType()); + } else if (stateFormula.hasBound()) { + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { + result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), storm::logic::OptimalityType::Maximize); + } else { + result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), storm::logic::OptimalityType::Minimize); + } } else { result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative); } @@ -186,6 +192,12 @@ namespace storm { std::unique_ptr result; if (stateFormula.hasOptimalityType()) { result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, stateFormula.getOptimalityType()); + } else if (stateFormula.hasBound()) { + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { + result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), storm::logic::OptimalityType::Maximize); + } else { + result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), storm::logic::OptimalityType::Minimize); + } } else { result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative); } @@ -212,6 +224,12 @@ namespace storm { std::unique_ptr result; if (stateFormula.hasOptimalityType()) { result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, stateFormula.getOptimalityType()); + } else if (stateFormula.hasBound()) { + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { + result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), storm::logic::OptimalityType::Maximize); + } else { + result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), storm::logic::OptimalityType::Minimize); + } } else { result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative); } @@ -227,13 +245,19 @@ 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 result; if (stateFormula.hasOptimalityType()) { - return this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, stateFormula.getOptimalityType()); + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, stateFormula.getOptimalityType()); + } else if (stateFormula.hasBound()) { + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), storm::logic::OptimalityType::Maximize); + } else { + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), storm::logic::OptimalityType::Minimize); + } } else { - return this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false); + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false); } - std::unique_ptr result; if (stateFormula.hasBound()) { return result->compareAgainstBound(stateFormula.getComparisonType(), stateFormula.getBound()); } else { diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 3aec3ef8d..ced5463a6 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -69,7 +69,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula()); std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); @@ -92,7 +92,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector()))); @@ -161,7 +161,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula()); std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); @@ -200,7 +200,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getStepBound()))); } @@ -220,7 +220,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getStepCount()))); } @@ -307,7 +307,7 @@ namespace storm { template std::unique_ptr SparseMdpPrctlModelChecker::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, 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 this->getModel()."); + STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic."); std::unique_ptr subResultPointer = this->check(rewardPathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 848abc452..630ce0bb2 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -10,6 +10,9 @@ namespace storm { namespace counterexamples { template class SMTMinimalCommandSetGenerator; + + template + class MILPMinimalLabelSetGenerator; } namespace modelchecker { @@ -22,7 +25,8 @@ namespace storm { class SparseMdpPrctlModelChecker : public SparsePropositionalModelChecker { public: friend class SparseMarkovAutomatonCslModelChecker; - friend class counterexamples::SMTMinimalCommandSetGenerator; + friend class storm::counterexamples::SMTMinimalCommandSetGenerator; + friend class storm::counterexamples::MILPMinimalLabelSetGenerator; explicit SparseMdpPrctlModelChecker(storm::models::Mdp const& model); explicit SparseMdpPrctlModelChecker(storm::models::Mdp const& model, std::shared_ptr> nondeterministicLinearEquationSolver); diff --git a/src/utility/cli.h b/src/utility/cli.h index e00c4f612..3e0c31b5a 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -321,6 +321,11 @@ namespace storm { } options.addConstantDefinitionsFromString(program, settings.getConstantDefinitionString()); + // Generate command labels if we are going to build a counterexample later. + if (storm::settings::counterexampleGeneratorSettings().isMinimalCommandSetGenerationSet()) { + options.buildCommandLabels = true; + } + // Then, build the model from the symbolic description. result = storm::builder::ExplicitPrismModelBuilder::translateProgram(program, options); return result; From 5bbd85c379802a6e57119e732328ce24bd98b5b4 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 17 Mar 2015 10:39:13 +0100 Subject: [PATCH 5/9] Some bugfixes. Former-commit-id: 70dcc73e910c4039a872a64b33c494d1f66b507c --- src/modelchecker/AbstractModelChecker.cpp | 24 +++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp index 6b4ab55ab..ed51c16d2 100644 --- a/src/modelchecker/AbstractModelChecker.cpp +++ b/src/modelchecker/AbstractModelChecker.cpp @@ -161,10 +161,10 @@ namespace storm { if (stateFormula.hasOptimalityType()) { result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, stateFormula.getOptimalityType()); } else if (stateFormula.hasBound()) { - if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { - result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), storm::logic::OptimalityType::Maximize); + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) { + result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, storm::logic::OptimalityType::Maximize); } else { - result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), storm::logic::OptimalityType::Minimize); + result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, storm::logic::OptimalityType::Minimize); } } else { result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative); @@ -193,10 +193,10 @@ namespace storm { if (stateFormula.hasOptimalityType()) { result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, stateFormula.getOptimalityType()); } else if (stateFormula.hasBound()) { - if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { - result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), storm::logic::OptimalityType::Maximize); + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) { + result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, storm::logic::OptimalityType::Maximize); } else { - result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), storm::logic::OptimalityType::Minimize); + result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, storm::logic::OptimalityType::Minimize); } } else { result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative); @@ -225,10 +225,10 @@ namespace storm { if (stateFormula.hasOptimalityType()) { result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, stateFormula.getOptimalityType()); } else if (stateFormula.hasBound()) { - if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { - result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), storm::logic::OptimalityType::Maximize); + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) { + result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, storm::logic::OptimalityType::Maximize); } else { - result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), storm::logic::OptimalityType::Minimize); + result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, storm::logic::OptimalityType::Minimize); } } else { result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative); @@ -249,10 +249,10 @@ namespace storm { if (stateFormula.hasOptimalityType()) { result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, stateFormula.getOptimalityType()); } else if (stateFormula.hasBound()) { - if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || storm::logic::ComparisonType::LessEqual) { - result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), storm::logic::OptimalityType::Maximize); + if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) { + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, storm::logic::OptimalityType::Maximize); } else { - result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), storm::logic::OptimalityType::Minimize); + result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, storm::logic::OptimalityType::Minimize); } } else { result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false); From 96539f41a500f518ec32d5c64cb416fcbdacd87d Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 17 Mar 2015 12:40:59 +0100 Subject: [PATCH 6/9] Fixed simplification of division: division expressions must not be simplified, because it is not (yet) clear whether integer division or floating point division is to be used. Former-commit-id: 506798c1cdc40d8311525661cd090888f724558b --- src/builder/ExplicitPrismModelBuilder.cpp | 13 +++++++++++-- .../BinaryNumericalFunctionExpression.cpp | 7 ++++--- .../UnaryNumericalFunctionExpression.cpp | 2 +- src/storage/prism/BooleanVariable.cpp | 2 +- 4 files changed, 17 insertions(+), 7 deletions(-) diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp index 5ac83af33..082cdb457 100644 --- a/src/builder/ExplicitPrismModelBuilder.cpp +++ b/src/builder/ExplicitPrismModelBuilder.cpp @@ -127,6 +127,7 @@ namespace storm { // all expressions in the program so we can then evaluate them without having to store the values of the // constants in the state (i.e., valuation). preparedProgram = preparedProgram.substituteConstants(); + storm::prism::RewardModel rewardModel = storm::prism::RewardModel(); // Select the appropriate reward model. @@ -564,7 +565,11 @@ namespace storm { } for (auto const& stateProbabilityPair : choice) { - globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices; + if (discreteTimeModel) { + globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices; + } else { + globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second; + } // Now add all rewards that match this choice. for (auto const& transitionReward : transitionRewards) { @@ -580,7 +585,11 @@ namespace storm { } for (auto const& stateProbabilityPair : choice) { - globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices; + if (discreteTimeModel) { + globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices; + } else { + globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second; + } // Now add all rewards that match this choice. for (auto const& transitionReward : transitionRewards) { diff --git a/src/storage/expressions/BinaryNumericalFunctionExpression.cpp b/src/storage/expressions/BinaryNumericalFunctionExpression.cpp index be0f30dc0..4e95ea551 100644 --- a/src/storage/expressions/BinaryNumericalFunctionExpression.cpp +++ b/src/storage/expressions/BinaryNumericalFunctionExpression.cpp @@ -6,6 +6,7 @@ #include "src/storage/expressions/DoubleLiteralExpression.h" #include "src/utility/macros.h" #include "src/exceptions/InvalidTypeException.h" +#include "src/exceptions/InvalidStateException.h" namespace storm { namespace expressions { @@ -65,7 +66,7 @@ namespace storm { std::shared_ptr firstOperandSimplified = this->getFirstOperand()->simplify(); std::shared_ptr secondOperandSimplified = this->getSecondOperand()->simplify(); - if (firstOperandSimplified->isLiteral() && secondOperandSimplified->isLiteral()) { + if (firstOperandSimplified->isLiteral() && secondOperandSimplified->isLiteral() && this->getOperatorType() != OperatorType::Divide) { if (this->hasIntegerType()) { int_fast64_t firstOperandEvaluation = firstOperandSimplified->evaluateAsInt(); int_fast64_t secondOperandEvaluation = secondOperandSimplified->evaluateAsInt(); @@ -74,10 +75,10 @@ namespace storm { case OperatorType::Plus: newValue = firstOperandEvaluation + secondOperandEvaluation; break; case OperatorType::Minus: newValue = firstOperandEvaluation - secondOperandEvaluation; break; case OperatorType::Times: newValue = firstOperandEvaluation * secondOperandEvaluation; break; - case OperatorType::Divide: newValue = firstOperandEvaluation / secondOperandEvaluation; break; case OperatorType::Min: newValue = std::min(firstOperandEvaluation, secondOperandEvaluation); break; case OperatorType::Max: newValue = std::max(firstOperandEvaluation, secondOperandEvaluation); break; case OperatorType::Power: newValue = static_cast(std::pow(firstOperandEvaluation, secondOperandEvaluation)); break; + case OperatorType::Divide: STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Unable to simplify division."); break; } return std::shared_ptr(new IntegerLiteralExpression(this->getManager(), newValue)); } else if (this->hasRationalType()) { @@ -88,10 +89,10 @@ namespace storm { case OperatorType::Plus: newValue = firstOperandEvaluation + secondOperandEvaluation; break; case OperatorType::Minus: newValue = firstOperandEvaluation - secondOperandEvaluation; break; case OperatorType::Times: newValue = firstOperandEvaluation * secondOperandEvaluation; break; - case OperatorType::Divide: newValue = firstOperandEvaluation / secondOperandEvaluation; break; case OperatorType::Min: newValue = std::min(firstOperandEvaluation, secondOperandEvaluation); break; case OperatorType::Max: newValue = std::max(firstOperandEvaluation, secondOperandEvaluation); break; case OperatorType::Power: newValue = static_cast(std::pow(firstOperandEvaluation, secondOperandEvaluation)); break; + case OperatorType::Divide: STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Unable to simplify division."); break; } return std::shared_ptr(new DoubleLiteralExpression(this->getManager(), newValue)); } diff --git a/src/storage/expressions/UnaryNumericalFunctionExpression.cpp b/src/storage/expressions/UnaryNumericalFunctionExpression.cpp index 6e2a13ae0..edfb22702 100644 --- a/src/storage/expressions/UnaryNumericalFunctionExpression.cpp +++ b/src/storage/expressions/UnaryNumericalFunctionExpression.cpp @@ -76,7 +76,7 @@ namespace storm { case OperatorType::Floor: value = std::floor(boost::get(operandEvaluation)); break; case OperatorType::Ceil: value = std::ceil(boost::get(operandEvaluation)); break; } - return std::shared_ptr(new DoubleLiteralExpression(this->getManager(), boost::get(result))); + return std::shared_ptr(new DoubleLiteralExpression(this->getManager(), value)); } } diff --git a/src/storage/prism/BooleanVariable.cpp b/src/storage/prism/BooleanVariable.cpp index 35df52356..c38812ec6 100644 --- a/src/storage/prism/BooleanVariable.cpp +++ b/src/storage/prism/BooleanVariable.cpp @@ -11,7 +11,7 @@ namespace storm { } std::ostream& operator<<(std::ostream& stream, BooleanVariable const& variable) { - stream << variable.getName() << ": bool init" << variable.getInitialValueExpression() << ";"; + stream << variable.getName() << ": bool init " << variable.getInitialValueExpression() << ";"; return stream; } From 3434405cf43c9a8e69600ee755eb6997e10e3714 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 17 Mar 2015 18:13:30 +0100 Subject: [PATCH 7/9] Started working on CTMC mc. Former-commit-id: 7e38c0d7d33d816a8cc397deedb48bb458dfcf1c --- src/utility/numerical.h | 158 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 src/utility/numerical.h diff --git a/src/utility/numerical.h b/src/utility/numerical.h new file mode 100644 index 000000000..eb9aaf48d --- /dev/null +++ b/src/utility/numerical.h @@ -0,0 +1,158 @@ +#include +#include +#include + +#include "src/utility/macros.h" +#include "src/utility/ConstantsComparator.h" +#include "src/exceptions/InvalidArgumentException.h" + +namespace storm { + namespace utility { + namespace numerical { + template + std::tuple> getFoxGlynnCutoff(ValueType lambda, ) { + storm::utility::ConstantsComparator comparator; + + ValueType underflow, overflow, accuracy = 0; + + STORM_LOG_THROW(!comparator.isZero(lambda), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero."); + + // This code is more or less taken from PRISM. According to their implementation, for lambda smaller than + // 400, we compute the result using the naive method. + if (lambda < 400) { + ValueType eToPowerMinusLambda = std::exp(-lambda); + ValueType targetValue = (1 - (this->accuracy / 2.0)) / eToPowerMinusLambda; + std::vector weights; + + ValueType exactlyKEvents = 1; + ValueType atMostKEvents = exactlyKEvents; + weights.push_back(exactlyKEvents * eToPowerMinusLambda); + + uint_fast64_t currentK = 1; + do { + exactlyKEvents *= lambda / k; + atMostKEvents += exactlyKEvents; + weights.push_back(exactlyKEvents * eToPowerMinusLambda); + ++k; + } while (atMostKEvents < targetValue); + + return std::make_tuple(0, k - 1, 1.0, weights); + } else { + STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10."); + + ValueType factor = 1e+10; + ValueType m = std::floor(lambda); + + // Now start the Finder algorithm to find the truncation points. + { + // Factors used by the corollaries explained in Fox & Glynns paper. + // Square root of pi. + Type sqrtpi = 1.77245385090551602729; + + // Square root of 2. + Type sqrt2 = 1.41421356237309504880; + + // Square root of lambda. + Type sqrtLambda = std::sqrt(lambda); + + Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda. + Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda. + } + } + + // Use Fox Glynn algorithm for lambda>400. + if (accuracy < 1e-10) { + LOG4CPLUS_ERROR(logger, "Given Value accuracy must at least be 1e-10."); + throw storm::exceptions::InvalidArgumentException("Error while computing FoxGlynn values. accuracy < 1e-10."); + } + // Factor from Fox&Glynns paper. The paper does not explain where it comes from. + Type factor = 1e+10; + + // Run the FINDER algorithm. + Type m = floor(lambda); + { + // Factores used by the corollaries explained in Fox&Glynns paper. + Type sqrtpi = 1.77245385090551602729; // Squareroot of PI. + Type sqrt2 = 1.41421356237309504880; // Squareroot of 2. + Type sqrtLambda = sqrt(lambda); // Sqareroot of Lambda. + Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda. + Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda. + + // Use Corollary 1 from the paper to find the right truncation point. + Type k = 4; + res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi); + /* Normally: Iterate between the bounds stated above to find the right truncationpoint. + * This is a modification to the Fox-Glynn paper. The search for the right truncationpoint is only + * terminated by the error condition and not by the upper bound. Thus the calculation can be more accurate. + */ + while (res > accuracy / 2.0) + { + k++; + Type dkl = 1.0 / (1 - exp(-(2.0 / 9.0)*(k*sqrt2*sqrtLambda + 1.5))); // d(k,Lambda) from the paper. + Type res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 1. + + } + + // Left hand side of the equation in Corollary 1. + this->rightTrunk = (int)ceil(m + k*sqrt2*sqrtLambda + 1.5); + + // Use Corollary 2 to find left truncation point. + Type res; + k = 4; + + do + { + res = b_Lambda*exp(-k*-k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 2 + k++; + } while (res > accuracy / 2.0); + + this->leftTrunk = (int)(m - k*sqrtLambda - 1.5); // Left hand side of the equation in Corollary 2. + + // Check for underflow. + if ((int)(m - k*sqrtLambda - 1.5) < 0.0) { + LOG4CPLUS_ERROR(logger, "Underflow while computing left truncation point."); + throw storm::exceptions::OutOfRangeException("Error while computing FoxGlynn values. Underflow of left Truncation point."); + } + + } + //End of FINDER algorithm. + + // Use WEIGHTER algorithm to determine weights. + // Down from m + for (Type j = m; j > this->leftTrunk; j--) + this->weights.at(j - 1 - this->leftTrunk) = (j / lambda)*this->weights.at(j - this->leftTrunk); + // Up from m + for (Type j = m; j < this->rightTrunk; j++) + this->weights.at(j + 1 - this->leftTrunk) = (lambda / (j + 1))*this->weights.at(j - this->leftTrunk); + + // Compute total weight. + // Add up weights from smallest to largest to avoid roundoff errors. + this->totalWeight = 0.0; + Type s = this->leftTrunk; + Type t = this->rightTrunk; + while (s < t) + { + if (this->weights.at(s - this->leftTrunk) <= this->weights.at(t - this->leftTrunk)) + { + this->totalWeight += this->weights.at(s - this->leftTrunk); + s++; + } + else + { + this->totalWeight += this->weights.at(t - this->leftTrunk); + t--; + } + } + + this->totalWeight += this->weights.at(s - this->leftTrunk); + + LOG4CPLUS_INFO(logger, "Left truncationpoint: " << this->leftTrunk << "."); + LOG4CPLUS_INFO(logger, "Right truncationpoint: " << this->rightTrunk << "."); + LOG4CPLUS_INFO(logger, "Total Weight:" << this->totalWeight << "."); + LOG4CPLUS_INFO(logger, "10. Weight: " << this->weights.at(9) << "."); + } + + } + } + } +} \ No newline at end of file From cde9786dfa0b76a47592cac64ed15b43ed72e5d1 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 17 Mar 2015 22:11:17 +0100 Subject: [PATCH 8/9] Made Fox-Glynn (hopefully) work. Former-commit-id: b07c88d122fd28642aaa24f0eb421849b9961224 --- src/utility/cli.h | 13 ++++ src/utility/numerical.h | 159 +++++++++++++++++----------------------- 2 files changed, 79 insertions(+), 93 deletions(-) diff --git a/src/utility/cli.h b/src/utility/cli.h index 3e0c31b5a..c493c5acd 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -75,6 +75,9 @@ log4cplus::Logger printer; #include "src/exceptions/InvalidSettingsException.h" #include "src/exceptions/InvalidTypeException.h" +// FIXME: remove this +#include "src/utility/numerical.h" + namespace storm { namespace utility { namespace cli { @@ -507,6 +510,16 @@ namespace storm { initializeFileLogging(); } + auto fgresult = storm::utility::numerical::getFoxGlynnCutoff(500, 1e-300, 1e+300, 1e-6); + + std::cout << "left: " << std::get<0>(fgresult) << " and right: " << std::get<1>(fgresult) << std::endl; + std::cout << "weight: " << std::get<2>(fgresult) << std::endl; + int pos = 0; + for (auto const& element : std::get<3>(fgresult)) { + std::cout << "elem[" << pos << "]: " << element << std::endl; + ++pos; + } + storm::settings::modules::GeneralSettings const& settings = storm::settings::generalSettings(); // If we have to build the model from a symbolic representation, we need to parse the representation first. diff --git a/src/utility/numerical.h b/src/utility/numerical.h index eb9aaf48d..d14036859 100644 --- a/src/utility/numerical.h +++ b/src/utility/numerical.h @@ -3,32 +3,29 @@ #include #include "src/utility/macros.h" -#include "src/utility/ConstantsComparator.h" +#include "src/utility/constants.h" #include "src/exceptions/InvalidArgumentException.h" namespace storm { namespace utility { namespace numerical { template - std::tuple> getFoxGlynnCutoff(ValueType lambda, ) { + std::tuple> getFoxGlynnCutoff(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) { storm::utility::ConstantsComparator comparator; - - ValueType underflow, overflow, accuracy = 0; - STORM_LOG_THROW(!comparator.isZero(lambda), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero."); - // This code is more or less taken from PRISM. According to their implementation, for lambda smaller than - // 400, we compute the result using the naive method. - if (lambda < 400) { + // This code is a modified version of the one in PRISM. According to their implementation, for lambda + // smaller than 400, we compute the result using the naive method. + if (lambda < 25) { ValueType eToPowerMinusLambda = std::exp(-lambda); - ValueType targetValue = (1 - (this->accuracy / 2.0)) / eToPowerMinusLambda; + ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda; std::vector weights; ValueType exactlyKEvents = 1; ValueType atMostKEvents = exactlyKEvents; weights.push_back(exactlyKEvents * eToPowerMinusLambda); - uint_fast64_t currentK = 1; + uint_fast64_t k = 1; do { exactlyKEvents *= lambda / k; atMostKEvents += exactlyKEvents; @@ -40,118 +37,94 @@ namespace storm { } else { STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10."); + // Factor from Fox&Glynn's paper. The paper does not explain where it comes from. ValueType factor = 1e+10; - ValueType m = std::floor(lambda); // Now start the Finder algorithm to find the truncation points. + ValueType m = std::floor(lambda); + uint_fast64_t leftTruncationPoint = 0, rightTruncationPoint = 0; { // Factors used by the corollaries explained in Fox & Glynns paper. // Square root of pi. - Type sqrtpi = 1.77245385090551602729; + ValueType sqrtpi = 1.77245385090551602729; // Square root of 2. - Type sqrt2 = 1.41421356237309504880; + ValueType sqrt2 = 1.41421356237309504880; - // Square root of lambda. - Type sqrtLambda = std::sqrt(lambda); - - Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda. - Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda. - } - } - - // Use Fox Glynn algorithm for lambda>400. - if (accuracy < 1e-10) { - LOG4CPLUS_ERROR(logger, "Given Value accuracy must at least be 1e-10."); - throw storm::exceptions::InvalidArgumentException("Error while computing FoxGlynn values. accuracy < 1e-10."); - } - // Factor from Fox&Glynns paper. The paper does not explain where it comes from. - Type factor = 1e+10; - - // Run the FINDER algorithm. - Type m = floor(lambda); - { - // Factores used by the corollaries explained in Fox&Glynns paper. - Type sqrtpi = 1.77245385090551602729; // Squareroot of PI. - Type sqrt2 = 1.41421356237309504880; // Squareroot of 2. - Type sqrtLambda = sqrt(lambda); // Sqareroot of Lambda. - Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda. - Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda. + // Set up a_\lambda, b_\lambda, and the square root of lambda. + ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0; + if (m < 400) { + sqrtLambda = std::sqrt(400.0); + aLambda = (1.0 + 1.0 / 400.0) * exp(0.0625) * sqrt2; + bLambda = (1.0 + 1.0 / 400.0) * exp(0.125 / 400.0); + } else { + sqrtLambda = std::sqrt(lambda); + aLambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; + bLambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); + } // Use Corollary 1 from the paper to find the right truncation point. - Type k = 4; - res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi); - /* Normally: Iterate between the bounds stated above to find the right truncationpoint. - * This is a modification to the Fox-Glynn paper. The search for the right truncationpoint is only - * terminated by the error condition and not by the upper bound. Thus the calculation can be more accurate. - */ - while (res > accuracy / 2.0) - { - k++; - Type dkl = 1.0 / (1 - exp(-(2.0 / 9.0)*(k*sqrt2*sqrtLambda + 1.5))); // d(k,Lambda) from the paper. - Type res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 1. + uint_fast64_t k = 3; + + ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5))); + + // According to David Jansen the upper bound can be ignored to achieve more accurate results. + // Right hand side of the equation in Corollary 1. + while ((accuracy / 2.0) < (aLambda * dkl * std::exp(-(k*k / 2.0)) / (k * sqrt2 * sqrtpi))) { + ++k; + // d(k,Lambda) from the paper. + dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5))); } // Left hand side of the equation in Corollary 1. - this->rightTrunk = (int)ceil(m + k*sqrt2*sqrtLambda + 1.5); + rightTruncationPoint = static_cast(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5))); // Use Corollary 2 to find left truncation point. - Type res; - k = 4; + k = 3; - do - { - res = b_Lambda*exp(-k*-k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 2 - k++; - } while (res > accuracy / 2.0); + // Right hand side of the equation in Corollary 2. + while ((accuracy / 2.0) < ((bLambda * exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) { + ++k; + } - this->leftTrunk = (int)(m - k*sqrtLambda - 1.5); // Left hand side of the equation in Corollary 2. + // Left hand side of the equation in Corollary 2. + leftTruncationPoint = std::max(static_cast(std::trunc(m - k * sqrtLambda - 1.5)), uint_fast64_t(0)); // Check for underflow. - if ((int)(m - k*sqrtLambda - 1.5) < 0.0) { - LOG4CPLUS_ERROR(logger, "Underflow while computing left truncation point."); - throw storm::exceptions::OutOfRangeException("Error while computing FoxGlynn values. Underflow of left Truncation point."); - } - + STORM_LOG_THROW(static_cast(std::trunc((m - k * sqrtLambda - 1.5))) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point."); } - //End of FINDER algorithm. - // Use WEIGHTER algorithm to determine weights. - // Down from m - for (Type j = m; j > this->leftTrunk; j--) - this->weights.at(j - 1 - this->leftTrunk) = (j / lambda)*this->weights.at(j - this->leftTrunk); - // Up from m - for (Type j = m; j < this->rightTrunk; j++) - this->weights.at(j + 1 - this->leftTrunk) = (lambda / (j + 1))*this->weights.at(j - this->leftTrunk); + std::vector weights(rightTruncationPoint - leftTruncationPoint + 1); + weights[m - leftTruncationPoint] = overflow / (factor * (rightTruncationPoint - leftTruncationPoint)); + for (uint_fast64_t j = m; j > leftTruncationPoint; --j) { + weights[j - 1 - leftTruncationPoint] = (j / lambda) * weights[j - leftTruncationPoint]; + } - // Compute total weight. - // Add up weights from smallest to largest to avoid roundoff errors. - this->totalWeight = 0.0; - Type s = this->leftTrunk; - Type t = this->rightTrunk; - while (s < t) - { - if (this->weights.at(s - this->leftTrunk) <= this->weights.at(t - this->leftTrunk)) - { - this->totalWeight += this->weights.at(s - this->leftTrunk); - s++; - } - else - { - this->totalWeight += this->weights.at(t - this->leftTrunk); - t--; + for (uint_fast64_t j = m; j < rightTruncationPoint; ++j) { + weights[j + 1 - leftTruncationPoint] = (lambda / (j + 1)) * weights[j - leftTruncationPoint]; + } + + + // Compute the total weight and start with smallest to avoid roundoff errors. + ValueType totalWeight = storm::utility::zero(); + + uint_fast64_t s = leftTruncationPoint; + uint_fast64_t t = rightTruncationPoint; + while (s < t) { + if (weights[s - leftTruncationPoint] <= weights[t - leftTruncationPoint]) { + totalWeight += weights[s - leftTruncationPoint]; + ++s; + } else { + totalWeight += weights[t - leftTruncationPoint]; + --t; } } - this->totalWeight += this->weights.at(s - this->leftTrunk); + totalWeight += weights[s - leftTruncationPoint]; - LOG4CPLUS_INFO(logger, "Left truncationpoint: " << this->leftTrunk << "."); - LOG4CPLUS_INFO(logger, "Right truncationpoint: " << this->rightTrunk << "."); - LOG4CPLUS_INFO(logger, "Total Weight:" << this->totalWeight << "."); - LOG4CPLUS_INFO(logger, "10. Weight: " << this->weights.at(9) << "."); + return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights); } - } } } From 9d4ded66b2583985f52844ac2b1a9b772f59348a Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 18 Mar 2015 11:52:45 +0100 Subject: [PATCH 9/9] Started implementing CTMC model checker. Former-commit-id: 8562e5e54c8f1698a800df8325fc57817e75c0d6 --- .../csl/SparseCtmcCslModelChecker.cpp | 88 +++++++++++++++++++ .../csl/SparseCtmcCslModelChecker.h | 44 ++++++++++ .../prctl/SparseDtmcPrctlModelChecker.cpp | 25 +++--- .../prctl/SparseDtmcPrctlModelChecker.h | 9 +- src/utility/cli.h | 5 -- 5 files changed, 150 insertions(+), 21 deletions(-) create mode 100644 src/modelchecker/csl/SparseCtmcCslModelChecker.cpp create mode 100644 src/modelchecker/csl/SparseCtmcCslModelChecker.h diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp new file mode 100644 index 000000000..56d48344d --- /dev/null +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -0,0 +1,88 @@ +#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" + +#include + +#include "src/utility/macros.h" +#include "src/utility/vector.h" +#include "src/utility/graph.h" +#include "src/utility/solver.h" + +#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" +#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" + +#include "src/exceptions/InvalidPropertyException.h" +#include "src/exceptions/NotImplementedException.h" + +namespace storm { + namespace modelchecker { + template + SparseCtmcCslModelChecker::SparseCtmcCslModelChecker(storm::models::Ctmc const& model) : SparsePropositionalModelChecker(model), linearEquationSolver(storm::utility::solver::getLinearEquationSolver()) { + // Intentionally left empty. + } + + template + SparseCtmcCslModelChecker::SparseCtmcCslModelChecker(storm::models::Ctmc const& model, std::unique_ptr>&& linearEquationSolver) : SparsePropositionalModelChecker(model), linearEquationSolver(std::move(linearEquationSolver)) { + // Intentionally left empty. + } + + template + bool SparseCtmcCslModelChecker::canHandle(storm::logic::Formula const& formula) const { + // FIXME: refine. + return formula.isCslStateFormula() || formula.isCslPathFormula() || formula.isRewardPathFormula(); + } + + template + std::unique_ptr SparseCtmcCslModelChecker::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { + STORM_LOG_THROW(pathFormula.isIntervalBounded(), storm::exceptions::InvalidPropertyException, "Cannot treat non-interval bounded until."); + + std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula()); + std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); + ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();; + ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); + std::pair const& intervalBounds = pathFormula.getIntervalBounds(); + std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), intervalBounds.first, intervalBounds.second))); + + return result; + } + + template + std::unique_ptr SparseCtmcCslModelChecker::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + std::vector result = SparseDtmcPrctlModelChecker::computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + } + + template + std::unique_ptr SparseCtmcCslModelChecker::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula()); + std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); + ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); + ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); + std::vector result = SparseDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + } + +// template +// std::unique_ptr SparseCtmcCslModelChecker::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional const& optimalityType) { +// +// } + + template + storm::models::Ctmc const& SparseCtmcCslModelChecker::getModel() const { + return this->template getModelAs>(); + } + + template + std::vector SparseCtmcCslModelChecker::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double lowerBound, double upperBound) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "not yet implemented, sorry"); + } + +// template +// std::vector SparseCtmcCslModelChecker::computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const { +// +// } + + } // namespace modelchecker +} // namespace storm \ No newline at end of file diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h new file mode 100644 index 000000000..7c5617a5e --- /dev/null +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -0,0 +1,44 @@ +#ifndef STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_ +#define STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_ + +#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" +#include "src/models/Ctmc.h" +#include "src/solver/LinearEquationSolver.h" + +namespace storm { + namespace modelchecker { + + template + class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker { + public: + explicit SparseCtmcCslModelChecker(storm::models::Ctmc const& model); + explicit SparseCtmcCslModelChecker(storm::models::Ctmc const& model, std::unique_ptr>&& linearEquationSolver); + + // The implemented methods of the AbstractModelChecker interface. + virtual bool canHandle(storm::logic::Formula const& formula) const override; + virtual std::unique_ptr computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; +// 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; + + protected: + storm::models::Ctmc const& getModel() const override; + + private: + // The methods that perform the actual checking. + std::vector computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double lowerBound, double upperBound) const; + std::vector computeUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const; +// 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; + + // An object that is used for solving linear equations and performing matrix-vector multiplication. + std::unique_ptr> linearEquationSolver; + }; + + } // namespace modelchecker +} // namespace storm + +#endif /* STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_ */ diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 9d12d6166..d3a01be83 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -5,6 +5,7 @@ #include "src/utility/macros.h" #include "src/utility/vector.h" #include "src/utility/graph.h" +#include "src/utility/solver.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" @@ -74,14 +75,13 @@ namespace storm { } template - std::vector SparseDtmcPrctlModelChecker::computeNextProbabilitiesHelper(storm::storage::BitVector const& nextStates) { + std::vector SparseDtmcPrctlModelChecker::computeNextProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver const& linearEquationSolver) { // Create the vector with which to multiply and initialize it correctly. - std::vector result(this->getModel().getNumberOfStates()); + std::vector result(transitionMatrix.getRowCount()); storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one()); // Perform one single matrix-vector multiplication. - STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available."); - this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), result); + linearEquationSolver.performMatrixVectorMultiplication(transitionMatrix, result); return result; } @@ -89,14 +89,14 @@ namespace storm { std::unique_ptr SparseDtmcPrctlModelChecker::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeNextProbabilitiesHelper(subResult.getTruthValuesVector()))); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver))); } template - std::vector SparseDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const { + std::vector SparseDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver const& linearEquationSolver) { // We need to identify the states which have to be taken out of the matrix, i.e. // all states that have probability 0 and 1 of satisfying the until-formula. - std::pair statesWithProbability01 = storm::utility::graph::performProb01(this->getModel(), phiStates, psiStates); + std::pair statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates); storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first); storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second); @@ -107,7 +107,7 @@ namespace storm { STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); // Create resulting vector. - std::vector result(this->getModel().getNumberOfStates()); + std::vector result(transitionMatrix.getRowCount()); // Check whether we need to compute exact probabilities for some states. if (qualitative) { @@ -118,7 +118,7 @@ namespace storm { // In this case we have have to compute the probabilities. // We can eliminate the rows and columns from the original transition probability matrix. - storm::storage::SparseMatrix submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, true); + storm::storage::SparseMatrix submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true); // Converting the matrix from the fixpoint notation to the form needed for the equation // system. That is, we go from x = A*x + b to (I-A)x = b. @@ -131,11 +131,10 @@ namespace storm { // Prepare the right-hand side of the equation system. For entry i this corresponds to // the accumulated probability of going from state i to some 'yes' state. - std::vector b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); + std::vector b = transitionMatrix.getConstrainedRowSumVector(maybeStates, statesWithProbability1); // Now solve the created system of linear equations. - STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available."); - this->linearEquationSolver->solveEquationSystem(submatrix, x, b); + linearEquationSolver.solveEquationSystem(submatrix, x, b); // Set values of resulting vector according to result. storm::utility::vector::setVectorValues(result, maybeStates, x); @@ -155,7 +154,7 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();; ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();; - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative))); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver))); } template diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index f7557d63b..9fe58e0f2 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -3,15 +3,18 @@ #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "src/models/Dtmc.h" -#include "src/utility/solver.h" #include "src/solver/LinearEquationSolver.h" namespace storm { namespace modelchecker { + // Forward-declare CTMC model checker so we can make it a friend. + template class SparseCtmcCslModelChecker; template class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker { public: + friend class SparseCtmcCslModelChecker; + explicit SparseDtmcPrctlModelChecker(storm::models::Dtmc const& model); explicit SparseDtmcPrctlModelChecker(storm::models::Dtmc const& model, std::unique_ptr>&& linearEquationSolver); @@ -30,8 +33,8 @@ namespace storm { private: // The methods that perform the actual checking. std::vector computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const; - std::vector computeNextProbabilitiesHelper(storm::storage::BitVector const& nextStates); - std::vector computeUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const; + static std::vector computeNextProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver const& linearEquationSolver); + static std::vector computeUntilProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver const& linearEquationSolver); 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; diff --git a/src/utility/cli.h b/src/utility/cli.h index c493c5acd..0486b13e3 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -75,9 +75,6 @@ log4cplus::Logger printer; #include "src/exceptions/InvalidSettingsException.h" #include "src/exceptions/InvalidTypeException.h" -// FIXME: remove this -#include "src/utility/numerical.h" - namespace storm { namespace utility { namespace cli { @@ -510,8 +507,6 @@ namespace storm { initializeFileLogging(); } - auto fgresult = storm::utility::numerical::getFoxGlynnCutoff(500, 1e-300, 1e+300, 1e-6); - std::cout << "left: " << std::get<0>(fgresult) << " and right: " << std::get<1>(fgresult) << std::endl; std::cout << "weight: " << std::get<2>(fgresult) << std::endl; int pos = 0;