From cde1c646d9359fdda973a2c1bea70743bb251888 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Mon, 11 Mar 2019 19:09:05 +0100
Subject: [PATCH 01/11] Started to implement the algorithm more close to the
one mentioned in the paper (in particular to make things more clean and to
allow more than 2 dimensions.
---
.../helper/rewardbounded/CostLimitClosure.cpp | 127 ++++++++++++
.../helper/rewardbounded/CostLimitClosure.h | 53 +++++
.../helper/rewardbounded/QuantileHelper.cpp | 184 +++++++++++++++++-
.../helper/rewardbounded/QuantileHelper.h | 7 +
4 files changed, 369 insertions(+), 2 deletions(-)
create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
new file mode 100644
index 000000000..622e674f1
--- /dev/null
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
@@ -0,0 +1,127 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
+#include "storm/utility/macros.h"
+#include "storm/exceptions/IllegalArgumentException.h"
+
+namespace storm {
+ namespace modelchecker {
+ namespace helper {
+ namespace rewardbounded {
+ CostLimit::CostLimit(uint64_t const &costLimit) : value(costLimit) {
+ // Intentionally left empty
+ }
+
+ bool CostLimit::isInfinity() const {
+ return value == std::numeric_limits<uint64_t>::max();
+ }
+
+ uint64_t const& CostLimit::get() const {
+ STORM_LOG_ASSERT(!isInfinity(), "Tried to get an infinite cost limit as int.");
+ return value;
+ }
+
+ uint64_t& CostLimit::get() {
+ STORM_LOG_ASSERT(!isInfinity(), "Tried to get an infinite cost limit as int.");
+ return value;
+ }
+
+ bool CostLimit::operator<(CostLimit const& other) const {
+ // Since infinity is represented by the max integer, we can compare this way.
+ return value < other.value;
+ }
+
+ bool CostLimit::operator==(CostLimit const& other) const {
+ return value == other.value;
+ }
+
+ CostLimit CostLimit::infinity() {
+ return CostLimit(std::numeric_limits<uint64_t>::max());
+ }
+
+ bool CostLimitClosure::CostLimitsCompare::operator()(storm::modelchecker::helper::rewardbounded::CostLimits const& lhs, storm::modelchecker::helper::rewardbounded::CostLimits const& rhs) const {
+ for (uint64_t i = 0; i < lhs.size(); ++i) {
+ if (lhs[i] < rhs[i]) {
+ return true;
+ } else if (rhs[i] < lhs[i]) {
+ return false;
+ }
+ }
+ return false;
+ }
+
+ CostLimitClosure::CostLimitClosure(storm::storage::BitVector const &downwardDimensions)
+ : downwardDimensions(downwardDimensions) {
+ // Intentionally left empty
+ }
+
+ bool CostLimitClosure::insert(CostLimits const& costLimits) {
+ // Iterate over all points in the generator and check whether they dominate the given point or vice versa
+ // TODO: make this more efficient by exploiting the order of the generator set.
+ std::vector<CostLimits> pointsToErase;
+ for (auto const& b : generator) {
+ if (dominates(b, costLimits)) {
+ // The given point is already contained in this closure.
+ // Since domination is transitive, this should not happen:
+ STORM_LOG_ASSERT(pointsToErase.empty(), "Inconsistent generator of CostLimitClosure.");
+ return false;
+ }
+ if (dominates(costLimits, b)) {
+ // b will be dominated by the new point so we erase it later.
+ // Note that b != newPoint holds if we reach this point
+ pointsToErase.push_back(b);
+ }
+ }
+ for (auto const& b : pointsToErase) {
+ generator.erase(b);
+ }
+ generator.insert(std::move(costLimits));
+ return true;
+ }
+
+ bool CostLimitClosure::contains(CostLimits const& costLimits) const {
+ // Iterate over all points in the generator and check whether they dominate the given point.
+ // TODO: make this more efficient by exploiting the order of the generator set.
+ for (auto const&b : generator) {
+ if (dominates(b,costLimits)) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ bool CostLimitClosure::dominates(CostLimits const& lhs, CostLimits const& rhs) const {
+ for (uint64_t i = 0; i < lhs.size(); ++i) {
+ if (downwardDimensions.get(i)) {
+ if (lhs[i] < rhs[i]) {
+ return false;
+ }
+ } else {
+ if (rhs[i] < lhs[i]) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+
+ std::vector<CostLimits> CostLimitClosure::getDominatingCostLimits(CostLimits const& costLimits) const {
+ std::vector<CostLimits> result;
+ for (auto const &b : generator) {
+ if (dominates(b, costLimits)) {
+ result.push_back(b);
+ }
+ }
+ return result;
+ }
+
+ typename CostLimitClosure::GeneratorType const &CostLimitClosure::getGenerator() const {
+ return generator;
+ }
+
+ uint64_t CostLimitClosure::dimension() const {
+ return downwardDimensions.size();
+ }
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
new file mode 100644
index 000000000..1da445d5f
--- /dev/null
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
@@ -0,0 +1,53 @@
+#pragma once
+
+#include <vector>
+#include <set>
+#include <string>
+
+#include "storm/storage/BitVector.h"
+
+namespace storm {
+ namespace modelchecker {
+ namespace helper {
+ namespace rewardbounded {
+
+ class CostLimit {
+ public:
+ CostLimit(uint64_t const& costLimit);
+ bool isInfinity() const;
+ uint64_t const& get() const;
+ uint64_t& get();
+ bool operator<(CostLimit const& other) const;
+ bool operator==(CostLimit const& other) const;
+ static CostLimit infinity();
+ private:
+ uint64_t value;
+ };
+ typedef std::vector<CostLimit> CostLimits;
+
+ class CostLimitClosure {
+ public:
+
+ struct CostLimitsCompare {
+ bool operator() (CostLimits const& lhs, CostLimits const& rhs) const;
+ };
+ typedef std::set<CostLimits, CostLimitsCompare> GeneratorType;
+
+ explicit CostLimitClosure(storm::storage::BitVector const& downwardDimensions);
+ bool insert(CostLimits const& costLimits);
+ bool contains(CostLimits const& costLimits) const;
+ bool dominates(CostLimits const& lhs, CostLimits const& rhs) const;
+ std::vector<CostLimits> getDominatingCostLimits(CostLimits const& costLimits) const;
+ GeneratorType const& getGenerator() const;
+ uint64_t dimension() const;
+
+ private:
+
+ /// The dimensions that are downwards closed, i.e., if x is in the closure, then also all y <= x
+ storm::storage::BitVector downwardDimensions;
+ GeneratorType generator;
+ };
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index 3a5ca0239..55c795aeb 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -39,6 +39,8 @@ namespace storm {
// Only > and <= are supported for upper bounds. This is to make sure that Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
// Only >= and < are supported for lower bounds. (EC construction..)
// TODO
+ // Fix precision increasing
+ // Multiple quantile formulas in the same file yield constants def clash
// Bounds are either constants or variables that are declared in the quantile formula.
// Prop op has optimality type
// No Prmin with lower cost bounds: Ec construction fails. In single obj we would do 1-Prmax[F "nonRewOrNonGoalEC"] but this invalidates other lower/upper cost bounds.
@@ -229,6 +231,184 @@ namespace storm {
return result;
}
+
+ template<typename ModelType>
+ std::pair<CostLimitClosure, std::vector<typename QuantileHelper<ModelType>::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions) const {
+ // Todo: ask the cache first
+
+ auto boundedUntilOp = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None));
+ std::set<storm::expressions::Variable> infinityVariables;
+ storm::storage::BitVector lowerBoundedDimensions(getDimension());
+ storm::storage::BitVector downwardClosedDimensions(getDimension());
+ for (auto const& d : getOpenDimensions()) {
+ if (consideredDimensions.get(d)) {
+ bool hasLowerCostBound = boundedUntilOp->getSubformula().asBoundedUntilFormula().hasLowerBound(d);
+ lowerBoundedDimensions.set(d, hasLowerCostBound);
+ bool hasLowerValueBound = storm::logic::isLowerBound(boundedUntilOp->getComparisonType());
+ downwardClosedDimensions.set(d, hasLowerCostBound == hasLowerValueBound);
+ } else {
+ infinityVariables.insert(getVariableForDimension(d));
+ }
+ }
+ downwardClosedDimensions = downwardClosedDimensions % consideredDimensions;
+ CostLimitClosure satCostLimits(downwardClosedDimensions), unsatCostLimits(~downwardClosedDimensions);
+
+ // Loop until the goal precision is reached.
+ while (true) {
+ // initialize Reward unfolding and data that will be needed for each epoch
+ MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, boundedUntilOp, infinityVariables);
+ if (computeQuantile(env, consideredDimensions, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) {
+ std::vector<ValueType> scalingFactors;
+ for (auto const& dim : consideredDimensions) {
+ scalingFactors.push_back(rewardUnfolding.getDimension(dim).scalingFactor);
+ }
+ return {satCostLimits, scalingFactors};
+ }
+ ++numPrecisionRefinements;
+ increasePrecision(env);
+ }
+ }
+
+
+ bool getNextCandidateCostLimit(CostLimit const& maxCostLimit, CostLimits& current) {
+ if (maxCostLimit.get() == 0) {
+ return false;
+ }
+ storm::storage::BitVector nonMaxEntries = storm::utility::vector::filter<CostLimit>(current, [&maxCostLimit] (CostLimit const& value) -> bool { return value < maxCostLimit; });
+ bool allZero = true;
+ for (auto const& entry : nonMaxEntries) {
+ if (current[entry].get() > 0) {
+ --current[entry].get();
+ allZero = false;
+ break;
+ } else {
+ current[entry] = CostLimit(maxCostLimit.get() - 1);
+ }
+ }
+ if (allZero) {
+ nonMaxEntries.increment();
+ if (nonMaxEntries.full()) {
+ return false;
+ }
+ current = CostLimits(current.size(), maxCostLimit);
+ storm::utility::vector::setVectorValues(current, nonMaxEntries, CostLimit(maxCostLimit.get() - 1));
+ }
+ return true;
+ }
+
+ bool unionContainsAllCostLimits(CostLimitClosure const& lhs, CostLimitClosure const& rhs) {
+ // TODO
+ assert(false);
+ }
+
+
+ bool translateEpochToCostLimits(EpochManager::Epoch const& epoch, EpochManager::Epoch const& startEpoch,storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, EpochManager const& epochManager, CostLimits& epochAsCostLimits) {
+ for (uint64_t dim = 0; dim < consideredDimensions.size(); ++dim) {
+ if (consideredDimensions.get(dim)) {
+ if (lowerBoundedDimensions.get(dim)) {
+ if (epochManager.isBottomDimension(epoch, dim)) {
+ epochAsCostLimits.push_back(CostLimit(0));
+ } else {
+ epochAsCostLimits.push_back(CostLimit(epochManager.getDimensionOfEpoch(epoch, dim) + 1));
+ }
+ } else {
+ if (epochManager.isBottomDimension(epoch, dim)) {
+ return false;
+ } else {
+ epochAsCostLimits.push_back(CostLimit(epochManager.getDimensionOfEpoch(epoch, dim)));
+ }
+ }
+ } else {
+ if (epochManager.isBottomDimension(epoch, dim)) {
+ if (!epochManager.isBottomDimension(startEpoch, dim)) {
+ return false;
+ }
+ } else if (epochManager.getDimensionOfEpoch(epoch, dim) != epochManager.getDimensionOfEpoch(startEpoch, dim)) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+
+ template<typename ModelType>
+ bool QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding) const {
+
+ auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+ auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+ std::vector<ValueType> x, b;
+ std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+ std::set<EpochManager::Epoch> checkedEpochs;
+
+
+ CostLimit currentMaxCostLimit(0);
+ for (CostLimit currentMaxCostLimit(0); !unionContainsAllCostLimits(satCostLimits, unsatCostLimits); ++currentMaxCostLimit.get()) {
+ CostLimits currentCandidate(satCostLimits.dimension(), currentMaxCostLimit);
+ do {
+ if (!satCostLimits.contains(currentCandidate) && !unsatCostLimits.contains(currentCandidate)) {
+ // Transform candidate cost limits to an appropriate start epoch
+ auto startEpoch = rewardUnfolding.getStartEpoch(true);
+ auto costLimitIt = currentCandidate.begin();
+ for (auto const& dim : consideredDimensions) {
+ if (lowerBoundedDimensions.get(dim)) {
+ if (costLimitIt->get() > 0) {
+ rewardUnfolding.getEpochManager().setDimensionOfEpoch(startEpoch, dim, costLimitIt->get() - 1);
+ } else {
+ rewardUnfolding.getEpochManager().setBottomDimension(startEpoch, dim);
+ }
+ } else {
+ rewardUnfolding.getEpochManager().setDimensionOfEpoch(startEpoch, dim, costLimitIt->get());
+ }
+ ++costLimitIt;
+ }
+ auto epochSequence = rewardUnfolding.getEpochComputationOrder(startEpoch);
+ for (auto const& epoch : epochSequence) {
+ if (checkedEpochs.count(epoch) == 0) {
+ checkedEpochs.insert(epoch);
+ auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
+ rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+ ++numCheckedEpochs;
+
+ CostLimits epochAsCostLimits;
+ if (translateEpochToCostLimits(epoch, startEpoch, consideredDimensions, lowerBoundedDimensions, rewardUnfolding.getEpochManager(), epochAsCostLimits)) {
+ ValueType currValue = rewardUnfolding.getInitialStateResult(epoch);
+ bool propertySatisfied;
+ if (env.solver().isForceSoundness()) {
+ auto lowerUpperValue = getLowerUpperBound(env, currValue);
+ propertySatisfied = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(lowerUpperValue.first);
+ if (propertySatisfied != quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(lowerUpperValue.second)) {
+ // unclear result due to insufficient precision.
+ return false;
+ }
+ } else {
+ propertySatisfied = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(currValue);
+ }
+
+ if (propertySatisfied) {
+ satCostLimits.insert(epochAsCostLimits);
+ } else {
+ unsatCostLimits.insert(epochAsCostLimits);
+ }
+ }
+ }
+ }
+ }
+ } while (getNextCandidateCostLimit(currentMaxCostLimit, currentCandidate));
+ }
+ return true;
+ }
+
+
+
+
+
+
+
+
+ ///////
+
+
template<typename ValueType>
void filterDominatedPoints(std::vector<std::vector<ValueType>>& points, std::vector<storm::solver::OptimizationDirection> const& dirs) {
std::vector<std::vector<ValueType>> result;
@@ -351,13 +531,13 @@ namespace storm {
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None)));
-
+
// initialize data that will be needed for each epoch
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
std::vector<ValueType> x, b;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
-
+
if (currentEpochValues[0] < 0 && currentEpochValues[1] < 0) {
// This case can only happen in these cases:
assert(lowerCostBounds[0]);
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index 13676a644..7b8b254cc 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -2,6 +2,8 @@
#include "storm/logic/QuantileFormula.h"
#include "boost/optional.hpp"
+#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
namespace storm {
class Environment;
@@ -23,6 +25,11 @@ namespace storm {
std::vector<std::vector<ValueType>> computeQuantile(Environment const& env) const;
private:
+
+ std::pair<CostLimitClosure, std::vector<ValueType>> computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions) const;
+ bool computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding);
+
+
std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env) const;
bool exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const;
From 8ae9a6f5d68479442a439c9689c1680d622dd6fe Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Tue, 12 Mar 2019 11:49:39 +0100
Subject: [PATCH 02/11] quantiles: Further improved the implementation as in
the paper
---
.../helper/rewardbounded/CostLimitClosure.cpp | 9 +-
.../helper/rewardbounded/CostLimitClosure.h | 1 +
.../MultiDimensionalRewardUnfolding.cpp | 53 +-
.../MultiDimensionalRewardUnfolding.h | 18 +-
.../helper/rewardbounded/ProductModel.cpp | 17 +-
.../prctl/helper/rewardbounded/ProductModel.h | 8 +-
.../helper/rewardbounded/QuantileHelper.cpp | 578 ++++--------------
.../helper/rewardbounded/QuantileHelper.h | 17 +-
8 files changed, 194 insertions(+), 507 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
index 622e674f1..f4fd31d85 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
@@ -88,6 +88,14 @@ namespace storm {
return false;
}
+ bool CostLimitClosure::containsUpwardClosure(CostLimits const& costLimits) const {
+ CostLimits infinityProjection(costLimits);
+ for (auto const& dim : downwardDimensions) {
+ infinityProjection[dim] = CostLimit::infinity();
+ }
+ return contains(infinityProjection);
+ }
+
bool CostLimitClosure::dominates(CostLimits const& lhs, CostLimits const& rhs) const {
for (uint64_t i = 0; i < lhs.size(); ++i) {
if (downwardDimensions.get(i)) {
@@ -103,7 +111,6 @@ namespace storm {
return true;
}
-
std::vector<CostLimits> CostLimitClosure::getDominatingCostLimits(CostLimits const& costLimits) const {
std::vector<CostLimits> result;
for (auto const &b : generator) {
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
index 1da445d5f..579a33b5b 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
@@ -36,6 +36,7 @@ namespace storm {
explicit CostLimitClosure(storm::storage::BitVector const& downwardDimensions);
bool insert(CostLimits const& costLimits);
bool contains(CostLimits const& costLimits) const;
+ bool containsUpwardClosure(CostLimits const& costLimits) const;
bool dominates(CostLimits const& lhs, CostLimits const& rhs) const;
std::vector<CostLimits> getDominatingCostLimits(CostLimits const& costLimits) const;
GeneratorType const& getGenerator() const;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
index 3b1f36791..5d82672eb 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
@@ -278,7 +278,8 @@ namespace storm {
STORM_LOG_THROW(SingleObjectiveMode, storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported in single objective mode."); // It most likely also works with multiple objectives with the same shape. However, we haven't checked this yet.
STORM_LOG_THROW(objectives.front().formula->isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for probability operator formulas");
auto const& probabilityOperatorFormula = objectives.front().formula->asProbabilityOperatorFormula();
- STORM_LOG_THROW(probabilityOperatorFormula.getSubformula().isBoundedUntilFormula() && probabilityOperatorFormula.hasOptimalityType() && storm::solver::maximize(probabilityOperatorFormula.getOptimalityType()), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for maximizing bounded until probabilities.");
+ STORM_LOG_THROW(probabilityOperatorFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for bounded until probabilities.");
+ STORM_LOG_THROW(!model.isNondeterministicModel() || (probabilityOperatorFormula.hasOptimalityType() && storm::solver::maximize(probabilityOperatorFormula.getOptimalityType())), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for maximizing bounded until probabilities.");
STORM_LOG_THROW(upperBoundedDimensions.empty() || !probabilityOperatorFormula.getSubformula().asBoundedUntilFormula().hasMultiDimensionalSubformulas(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported if the formula has either only lower bounds or if it has a single goal state."); // This would fail because the upper bounded dimension(s) might be satisfied already. One should erase epoch steps in the epoch model (after applying the goal-unfolding).
storm::storage::BitVector choicesWithoutUpperBoundedStep(model.getNumberOfChoices(), true);
@@ -373,7 +374,6 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
EpochModel<ValueType, SingleObjectiveMode>& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpoch(Epoch const& epoch) {
STORM_LOG_DEBUG("Setting model for epoch " << epochManager.toString(epoch));
- STORM_LOG_THROW(getProb1Objectives().empty(), storm::exceptions::InvalidOperationException, "The objective " << objectives[*getProb1Objectives().begin()].formula << " is satisfied already in the initial state(s). This special case can not be handled by the reward unfolding."); // This case should have been handled 'from the outside'
// Check if we need to update the current epoch class
if (!currentEpoch || !epochManager.compareEpochClass(epoch, currentEpoch.get())) {
@@ -669,6 +669,20 @@ namespace storm {
storm::utility::vector::addScaledVector(solution, solutionToAdd, scalingFactor);
}
+ template<typename ValueType, bool SingleObjectiveMode>
+ template<bool SO, typename std::enable_if<SO, int>::type>
+ void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const {
+ STORM_LOG_ASSERT(objIndex == 0, "Invalid objective index in single objective mode");
+ solution = value;
+ }
+
+ template<typename ValueType, bool SingleObjectiveMode>
+ template<bool SO, typename std::enable_if<!SO, int>::type>
+ void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const {
+ STORM_LOG_ASSERT(objIndex < solution.size(), "Invalid objective index " << objIndex << ".");
+ solution[objIndex] = value;
+ }
+
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<SO, int>::type>
std::string MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::solutionToString(SolutionType const& solution) const {
@@ -858,21 +872,39 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState) {
STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution at an unexisting product state.");
- STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at a state for which no solution was stored.");
- STORM_LOG_ASSERT(getProb1Objectives().empty(), "One of the objectives is already satisfied in the initial state. This special case is not handled."); // This case should have been handled 'from the outside'
+ STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at product state " << productState << " for which no solution was stored.");
return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]];
}
template<typename ValueType, bool SingleObjectiveMode>
- typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch) {
+ typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch) {
STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "The model has multiple initial states.");
- return getStateSolution(epoch, productModel->getInitialProductState(*model.getInitialStates().begin(), model.getInitialStates()));
+ return getInitialStateResult(epoch, *model.getInitialStates().begin());
}
template<typename ValueType, bool SingleObjectiveMode>
- typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) {
+ typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) {
STORM_LOG_ASSERT(model.getInitialStates().get(initialStateIndex), "The given model state is not an initial state.");
- return getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates()));
+
+ auto result = getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates(), epochManager.getEpochClass(epoch)));
+ for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+ if (productModel->getProb1InitialStates(objIndex) && productModel->getProb1InitialStates(objIndex)->get(initialStateIndex)) {
+ // Check whether the objective can actually hold in this epoch
+ bool objectiveHolds = true;
+ for (auto const& dim : objectiveDimensions[objIndex]) {
+ if (dimensions[dim].boundType == DimensionBoundType::LowerBound && !epochManager.isBottomDimension(epoch, dim)) {
+ objectiveHolds = false;
+ } else if (dimensions[dim].boundType == DimensionBoundType::UpperBound && epochManager.isBottomDimension(epoch, dim)) {
+ objectiveHolds = false;
+ }
+ STORM_LOG_ASSERT(dimensions[dim].boundType != DimensionBoundType::LowerBoundInfinity, "Unexpected bound type at this point.");
+ }
+ if (objectiveHolds) {
+ setSolutionEntry(result, objIndex, storm::utility::one<ValueType>());
+ }
+ }
+ }
+ return result;
}
template<typename ValueType, bool SingleObjectiveMode>
@@ -885,11 +917,6 @@ namespace storm {
return dimensions.at(dim);
}
- template<typename ValueType, bool SingleObjectiveMode>
- storm::storage::BitVector const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getProb1Objectives() const {
- return productModel->getProb1Objectives();
- }
-
template class MultiDimensionalRewardUnfolding<double, true>;
template class MultiDimensionalRewardUnfolding<double, false>;
template class MultiDimensionalRewardUnfolding<storm::RationalNumber, true>;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
index 510cdcfd5..3844c1af7 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
@@ -72,18 +72,12 @@ namespace storm {
boost::optional<ValueType> getLowerObjectiveBound(uint64_t objectiveIndex = 0);
void setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions);
- SolutionType const& getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
- SolutionType const& getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);
+ SolutionType getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
+ SolutionType getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);
EpochManager const& getEpochManager() const;
Dimension<ValueType> const& getDimension(uint64_t dim) const;
- /*!
- * Returns objectives that are always satisfied (i.e., have probability one) in all initial states.
- * These objectives can not be handled by this as they can not be translated into expected rewards.
- */
- storm::storage::BitVector const& getProb1Objectives() const;
-
private:
void setCurrentEpochClass(Epoch const& epoch);
@@ -105,6 +99,14 @@ namespace storm {
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const;
+
+ template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
+ void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const;
+ template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
+ void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const;
+
+
+
template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
std::string solutionToString(SolutionType const& solution) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
index aa387d084..ba86a26f9 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
@@ -21,7 +21,7 @@ namespace storm {
namespace rewardbounded {
template<typename ValueType>
- ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()), prob1Objectives(objectives.size(), false) {
+ ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()), prob1InitialStates(objectives.size(), boost::none) {
for (uint64_t dim = 0; dim < dimensions.size(); ++dim) {
if (!dimensions[dim].memoryLabel) {
@@ -167,10 +167,9 @@ namespace storm {
if (memStateBV.full()) {
storm::storage::BitVector initialTransitionStates = model.getInitialStates() & transitionStates;
// At this point we can check whether there is an initial state that already satisfies all subObjectives.
- // Such a situation can not be reduced (easily) to an expected reward computation.
- if (memStatePrimeBV.empty() && !initialTransitionStates.empty() && !objectiveContainsLowerBound && !initialTransitionStates.isDisjointFrom(constraintStates)) {
- STORM_LOG_THROW(model.getInitialStates() == initialTransitionStates, storm::exceptions::NotSupportedException, "The objective " << *objectives[objIndex].formula << " is already satisfied in an initial state. This special case is not supported.");
- prob1Objectives.set(objIndex, true);
+ // Such a situation can not be reduced (easily) to an expected reward computation and thus requires special treatment
+ if (memStatePrimeBV.empty() && !initialTransitionStates.empty()) {
+ prob1InitialStates[objIndex] = initialTransitionStates;
}
for (auto const& initState : initialTransitionStates) {
@@ -331,11 +330,11 @@ namespace storm {
}
template<typename ValueType>
- uint64_t ProductModel<ValueType>::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const {
+ uint64_t ProductModel<ValueType>::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates, EpochClass const& epochClass) const {
auto productInitStateIt = getProduct().getInitialStates().begin();
productInitStateIt += initialModelStates.getNumberOfSetBitsBeforeIndex(initialModelState);
STORM_LOG_ASSERT(getModelState(*productInitStateIt) == initialModelState, "Could not find the corresponding initial state in the product model.");
- return transformProductState(*productInitStateIt, epochManager.getEpochClass(epochManager.getZeroEpoch()), memoryStateManager.getInitialMemoryState());
+ return transformProductState(*productInitStateIt, epochClass, memoryStateManager.getInitialMemoryState());
}
template<typename ValueType>
@@ -649,8 +648,8 @@ namespace storm {
}
template<typename ValueType>
- storm::storage::BitVector const& ProductModel<ValueType>::getProb1Objectives() {
- return prob1Objectives;
+ boost::optional<storm::storage::BitVector> const& ProductModel<ValueType>::getProb1InitialStates(uint64_t objectiveIndex) const {
+ return prob1InitialStates[objectiveIndex];
}
template class ProductModel<double>;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
index 6ae408930..219d128bc 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
@@ -33,7 +33,7 @@ namespace storm {
bool productStateExists(uint64_t const& modelState, uint64_t const& memoryState) const;
uint64_t getProductState(uint64_t const& modelState, uint64_t const& memoryState) const;
- uint64_t getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const;
+ uint64_t getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates, EpochClass const& epochClass) const;
uint64_t getModelState(uint64_t const& productState) const;
MemoryState getMemoryState(uint64_t const& productState) const;
MemoryStateManager const& getMemoryStateManager() const;
@@ -47,8 +47,8 @@ namespace storm {
MemoryState transformMemoryState(MemoryState const& memoryState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
uint64_t transformProductState(uint64_t const& productState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
- // returns objectives that have probability one, already in the initial state.
- storm::storage::BitVector const& getProb1Objectives();
+ /// returns the initial states (with respect to the original model) that already satisfy the given objective with probability one, assuming that the cost bounds at the current epoch allow for the objective to be satisfied.
+ boost::optional<storm::storage::BitVector> const& getProb1InitialStates(uint64_t objectiveIndex) const;
private:
@@ -77,7 +77,7 @@ namespace storm {
std::vector<uint64_t> productToModelStateMap;
std::vector<MemoryState> productToMemoryStateMap;
std::vector<uint64_t> choiceToStateMap;
- storm::storage::BitVector prob1Objectives; /// Objectives that are already satisfied in the initial state
+ std::vector<boost::optional<storm::storage::BitVector>> prob1InitialStates; /// For each objective the set of initial states that already satisfy the objective
};
}
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index 55c795aeb..e78354752 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -29,27 +29,44 @@ namespace storm {
template<typename ModelType>
QuantileHelper<ModelType>::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) {
- // Intentionally left empty
+ // Do all kinds of sanity check.
+ std::set<storm::expressions::Variable> quantileVariables;
+ for (auto const& quantileVariable : quantileFormula.getBoundVariables()) {
+ STORM_LOG_THROW(quantileVariables.count(quantileVariable.second) == 0, storm::exceptions::NotSupportedException, "Quantile formula considers the same bound variable twice.");
+ quantileVariables.insert(quantileVariable.second);
+ }
STORM_LOG_THROW(quantileFormula.getSubformula().isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs probability operator inside. The formula " << quantileFormula << " is not supported.");
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
- STORM_LOG_THROW(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, storm::exceptions::NotSupportedException, "Probability operator inside quantile formula needs to have bound > or <=.");
+ STORM_LOG_THROW(probOpFormula.hasBound(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have a bound.");
+ STORM_LOG_THROW(!model.isNondeterministicModel() || probOpFormula.hasOptimalityType(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have an optimality type.");
+ STORM_LOG_WARN_COND(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, "Probability operator inside quantile formula needs to have bound > or <=. The specified comparison type might lead to non-termination."); // This has to do with letting bound variables approach infinity, e.g., Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
+ bool lowerBounded = storm::logic::isLowerBound(probOpFormula.getBound().comparisonType);
STORM_LOG_THROW(probOpFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs bounded until probability operator formula as subformula. The formula " << quantileFormula << " is not supported.");
- auto const& boundedUntilFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
+ auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
+ std::set<storm::expressions::Variable> boundVariables;
+ for (uint64_t dim = 0; dim < boundedUntilFormula.getDimension(); ++dim) {
+ storm::expressions::Expression boundExpression;
+ if (boundedUntilFormula.hasUpperBound(dim)) {
+ STORM_LOG_THROW(!boundedUntilFormula.hasLowerBound(dim), storm::exceptions::NotSupportedException, "Interval bounds are not supported within quantile formulas.");
+ STORM_LOG_THROW(!boundedUntilFormula.isUpperBoundStrict(dim), storm::exceptions::NotSupportedException, "Only non-strict upper reward bounds are supported for quantiles.");
+ boundExpression = boundedUntilFormula.getUpperBound(dim);
+ } else if (boundedUntilFormula.hasLowerBound(dim)) {
+ STORM_LOG_THROW(!boundedUntilFormula.isLowerBoundStrict(dim), storm::exceptions::NotSupportedException, "Only non-strict lower reward bounds are supported for quantiles.");
+ boundExpression = boundedUntilFormula.getLowerBound(dim);
+ }
+ if (boundExpression.isInitialized() && boundExpression.containsVariables()) {
+ STORM_LOG_THROW(boundExpression.isVariable(), storm::exceptions::NotSupportedException, "Non-trivial bound expressions such as '" << boundExpression << "' are not supported. Either specify a constant or a quantile variable.");
+ storm::expressions::Variable const& boundVariable = boundExpression.getBaseExpression().asVariableExpression().getVariable();
+ STORM_LOG_THROW(boundVariables.count(boundVariable) == 0, storm::exceptions::NotSupportedException, "Variable " << boundExpression << " occurs at multiple reward bounds.");
+ boundVariables.insert(boundVariable);
+ STORM_LOG_THROW(quantileVariables.count(boundVariable) == 1, storm::exceptions::NotSupportedException, "The formula contains undefined constant '" << boundExpression << "'.");
+ }
+ }
- // Only > and <= are supported for upper bounds. This is to make sure that Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
- // Only >= and < are supported for lower bounds. (EC construction..)
// TODO
// Fix precision increasing
// Multiple quantile formulas in the same file yield constants def clash
- // Bounds are either constants or variables that are declared in the quantile formula.
- // Prop op has optimality type
- // No Prmin with lower cost bounds: Ec construction fails. In single obj we would do 1-Prmax[F "nonRewOrNonGoalEC"] but this invalidates other lower/upper cost bounds.
- /* Todo: Current assumptions:
- * Subformula is always prob operator with bounded until
- * Each bound variable occurs at most once (change this?)
- * cost bounds are assumed to be always non-strict (the epochs returned by the reward unfolding assume strict lower and non-strict upper cost bounds, though...)
- * 'reasonable' quantile (e.g. not quantile(max B, Pmax>0.5 [F <=B G]) (we could also filter out these things early on)
- */
+ // ignore optimization direction for quantile variables
}
enum class BoundTransformation {
@@ -58,7 +75,7 @@ namespace storm {
GreaterEqualZero,
LessEqualZero
};
- std::shared_ptr<storm::logic::ProbabilityOperatorFormula> transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector<BoundTransformation> const& transformations) {
+ std::shared_ptr<storm::logic::ProbabilityOperatorFormula> transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector<BoundTransformation> const& transformations, bool complementQuery = false) {
auto const& origBoundedUntil = boundedUntilOperator.getSubformula().asBoundedUntilFormula();
STORM_LOG_ASSERT(transformations.size() == origBoundedUntil.getDimension(), "Tried to replace the bound of a dimension that is higher than the number of dimensions of the formula.");
std::vector<std::shared_ptr<storm::logic::Formula const>> leftSubformulas, rightSubformulas;
@@ -107,7 +124,11 @@ namespace storm {
} else {
newBoundedUntil = std::make_shared<storm::logic::BoundedUntilFormula>(origBoundedUntil.getLeftSubformula().asSharedPointer(), origBoundedUntil.getRightSubformula().asSharedPointer(), lowerBounds, upperBounds, timeBoundReferences);
}
- return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, boundedUntilOperator.getOperatorInformation());
+ storm::logic::OperatorInformation newOpInfo(boundedUntilOperator.getOperatorInformation().optimalityType, boundedUntilOperator.getBound());
+ if (complementQuery) {
+ newOpInfo.bound->comparisonType = storm::logic::invert(newOpInfo.bound->comparisonType);
+ }
+ return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, newOpInfo);
}
/// Increases the precision of solver results
@@ -174,55 +195,26 @@ namespace storm {
}
template<typename ModelType>
- storm::storage::BitVector QuantileHelper<ModelType>::getDimensionsForVariable(storm::expressions::Variable const& var) const {
- auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
- storm::storage::BitVector result(boundedUntil.getDimension(), false);
- for (uint64_t dim = 0; dim < boundedUntil.getDimension(); ++dim) {
- if (boundedUntil.hasLowerBound(dim) && boundedUntil.getLowerBound(dim).isVariable() && boundedUntil.getLowerBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) {
- result.set(dim, true);
- }
- if (boundedUntil.hasUpperBound(dim) && boundedUntil.getUpperBound(dim).isVariable() && boundedUntil.getUpperBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) {
- result.set(dim, true);
- }
- }
- return result;
- }
-
- template<typename ModelType>
- std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) const {
+ std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) {
+ numCheckedEpochs = 0;
+ numPrecisionRefinements = 0;
+ cachedSubQueryResults.clear();
+
std::vector<std::vector<ValueType>> result;
Environment envCpy = env; // It might be necessary to increase the precision during the computation
- numCheckedEpochs = 0; numPrecisionRefinements = 0;
- if (getOpenDimensions().getNumberOfSetBits() == 1) {
- uint64_t dimension = *getOpenDimensions().begin();
-
- bool zeroSatisfiesFormula = checkLimitValue(envCpy, storm::storage::BitVector(getDimension(), false));
- bool infSatisfiesFormula = checkLimitValue(envCpy, getOpenDimensions());
- if (zeroSatisfiesFormula != infSatisfiesFormula) {
- while (true) {
- auto quantileRes = computeQuantileForDimension(envCpy, dimension);
- if (quantileRes) {
- result = {{storm::utility::convertNumber<ValueType>(quantileRes->first) * quantileRes->second}};
- break;
- }
- increasePrecision(envCpy);
- ++numPrecisionRefinements;
- }
- } else if (zeroSatisfiesFormula) {
- // thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula
- if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) {
- result = {{ storm::utility::zero<ValueType>() }};
+ // Call the internal recursive function
+ auto internalResult = computeQuantile(envCpy, getOpenDimensions(), false);
+
+ // Translate the result by applying the scaling factors.
+ for (auto const& costLimits : internalResult.first.getGenerator()) {
+ std::vector<ValueType> resultPoint(costLimits.size());
+ storm::utility::vector::applyPointwise<CostLimit, ValueType, ValueType>(costLimits, internalResult.second, resultPoint, [](CostLimit const& costLimit, ValueType const& factor) -> ValueType {
+ if (costLimit.isInfinity()) {
+ return storm::utility::infinity<ValueType>();
} else {
- result = {{ storm::utility::infinity<ValueType>()}};
- }
- } else {
- // i.e., no bound value satisfies the formula
- result = {{}};
- }
- } else if (getOpenDimensions().getNumberOfSetBits() == 2) {
- result = computeTwoDimensionalQuantile(envCpy);
- } else {
- STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The quantile formula considers " << getOpenDimensions().getNumberOfSetBits() << " open dimensions. Only one- or two-dimensional quantiles are supported.");
+ return storm::utility::convertNumber<ValueType>(costLimit.get()) * factor;
+ }});
+ result.push_back(resultPoint);
}
if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl;
@@ -231,20 +223,26 @@ namespace storm {
return result;
}
-
template<typename ModelType>
- std::pair<CostLimitClosure, std::vector<typename QuantileHelper<ModelType>::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions) const {
- // Todo: ask the cache first
-
- auto boundedUntilOp = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None));
+ std::pair<CostLimitClosure, std::vector<typename QuantileHelper<ModelType>::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, bool complementaryQuery) {
+ STORM_LOG_ASSERT(consideredDimensions.isSubsetOf(getOpenDimensions()), "Considered dimensions for a quantile query should be a subset of the set of dimensions without a fixed bound.");
+
+ storm::storage::BitVector cacheKey = consideredDimensions;
+ cacheKey.resize(cacheKey.size() + 1, complementaryQuery);
+ auto cacheIt = cachedSubQueryResults.find(cacheKey);
+ if (cacheIt != cachedSubQueryResults.end()) {
+ return cacheIt->second;
+ }
+
+ auto boundedUntilOp = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None), complementaryQuery);
std::set<storm::expressions::Variable> infinityVariables;
storm::storage::BitVector lowerBoundedDimensions(getDimension());
storm::storage::BitVector downwardClosedDimensions(getDimension());
+ bool hasLowerValueBound = storm::logic::isLowerBound(boundedUntilOp->getComparisonType());
for (auto const& d : getOpenDimensions()) {
if (consideredDimensions.get(d)) {
bool hasLowerCostBound = boundedUntilOp->getSubformula().asBoundedUntilFormula().hasLowerBound(d);
lowerBoundedDimensions.set(d, hasLowerCostBound);
- bool hasLowerValueBound = storm::logic::isLowerBound(boundedUntilOp->getComparisonType());
downwardClosedDimensions.set(d, hasLowerCostBound == hasLowerValueBound);
} else {
infinityVariables.insert(getVariableForDimension(d));
@@ -253,23 +251,57 @@ namespace storm {
downwardClosedDimensions = downwardClosedDimensions % consideredDimensions;
CostLimitClosure satCostLimits(downwardClosedDimensions), unsatCostLimits(~downwardClosedDimensions);
+ // Initialize the (un)sat cost limits to guarantee termination
+ bool onlyUpperCostBounds = lowerBoundedDimensions.empty();
+ bool onlyLowerCostBounds = lowerBoundedDimensions == consideredDimensions;
+ if (onlyUpperCostBounds || onlyLowerCostBounds) {
+ for (auto const& k : consideredDimensions) {
+ storm::storage::BitVector subQueryDimensions = consideredDimensions;
+ subQueryDimensions.set(k, false);
+ bool subQueryComplement = complementaryQuery != ((onlyUpperCostBounds && hasLowerValueBound) || (onlyLowerCostBounds && !hasLowerValueBound));
+ auto subQueryResult = computeQuantile(env, subQueryDimensions, subQueryComplement);
+ for (auto const& subQueryCostLimit : subQueryResult.first.getGenerator()) {
+ CostLimits initPoint;
+ uint64_t i = 0;
+ for (auto const& dim : consideredDimensions) {
+ if (dim == k) {
+ initPoint.push_back(CostLimit::infinity());
+ } else {
+ initPoint.push_back(subQueryCostLimit[i]);
+ ++i;
+ }
+ }
+ if (subQueryComplement) {
+ unsatCostLimits.insert(initPoint);
+ } else {
+ satCostLimits.insert(initPoint);
+ }
+ }
+ }
+ } else {
+ STORM_LOG_WARN("Quantile formula considers mixtures of upper and lower reward-bounds. Termination is not guaranteed.");
+ }
+
// Loop until the goal precision is reached.
+ STORM_LOG_DEBUG("Computing quantile for dimensions: " << consideredDimensions);
while (true) {
- // initialize Reward unfolding and data that will be needed for each epoch
+ // initialize reward unfolding and data that will be needed for each epoch
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, boundedUntilOp, infinityVariables);
- if (computeQuantile(env, consideredDimensions, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) {
+ if (computeQuantile(env, consideredDimensions, *boundedUntilOp, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) {
std::vector<ValueType> scalingFactors;
for (auto const& dim : consideredDimensions) {
scalingFactors.push_back(rewardUnfolding.getDimension(dim).scalingFactor);
}
- return {satCostLimits, scalingFactors};
+ std::pair<CostLimitClosure, std::vector<ValueType>> result(satCostLimits, scalingFactors);
+ cachedSubQueryResults.emplace(cacheKey, result);
+ return result;
}
+ STORM_LOG_WARN("Restarting quantile computation due to insufficient precision.");
++numPrecisionRefinements;
increasePrecision(env);
}
}
-
bool getNextCandidateCostLimit(CostLimit const& maxCostLimit, CostLimits& current) {
if (maxCostLimit.get() == 0) {
return false;
@@ -296,12 +328,6 @@ namespace storm {
return true;
}
- bool unionContainsAllCostLimits(CostLimitClosure const& lhs, CostLimitClosure const& rhs) {
- // TODO
- assert(false);
- }
-
-
bool translateEpochToCostLimits(EpochManager::Epoch const& epoch, EpochManager::Epoch const& startEpoch,storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, EpochManager const& epochManager, CostLimits& epochAsCostLimits) {
for (uint64_t dim = 0; dim < consideredDimensions.size(); ++dim) {
if (consideredDimensions.get(dim)) {
@@ -331,9 +357,8 @@ namespace storm {
return true;
}
-
template<typename ModelType>
- bool QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding) const {
+ bool QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding) {
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
@@ -341,12 +366,14 @@ namespace storm {
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
std::set<EpochManager::Epoch> checkedEpochs;
-
- CostLimit currentMaxCostLimit(0);
- for (CostLimit currentMaxCostLimit(0); !unionContainsAllCostLimits(satCostLimits, unsatCostLimits); ++currentMaxCostLimit.get()) {
+ bool progress = true;
+ for (CostLimit currentMaxCostLimit(0); progress; ++currentMaxCostLimit.get()) {
CostLimits currentCandidate(satCostLimits.dimension(), currentMaxCostLimit);
+ // We can only stop the exploration, if the upward closure of the point in the 'top right corner' is contained in the (un)satCostlLimits.
+ progress = !satCostLimits.containsUpwardClosure(currentCandidate) && !unsatCostLimits.containsUpwardClosure(currentCandidate);
do {
if (!satCostLimits.contains(currentCandidate) && !unsatCostLimits.contains(currentCandidate)) {
+ progress = true;
// Transform candidate cost limits to an appropriate start epoch
auto startEpoch = rewardUnfolding.getStartEpoch(true);
auto costLimitIt = currentCandidate.begin();
@@ -367,7 +394,7 @@ namespace storm {
if (checkedEpochs.count(epoch) == 0) {
checkedEpochs.insert(epoch);
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+ rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env,boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
++numCheckedEpochs;
CostLimits epochAsCostLimits;
@@ -376,15 +403,14 @@ namespace storm {
bool propertySatisfied;
if (env.solver().isForceSoundness()) {
auto lowerUpperValue = getLowerUpperBound(env, currValue);
- propertySatisfied = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(lowerUpperValue.first);
- if (propertySatisfied != quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(lowerUpperValue.second)) {
+ propertySatisfied = boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.first);
+ if (propertySatisfied != boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.second)) {
// unclear result due to insufficient precision.
return false;
}
} else {
- propertySatisfied = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound().isSatisfied(currValue);
+ propertySatisfied = boundedUntilOperator.getBound().isSatisfied(currValue);
}
-
if (propertySatisfied) {
satCostLimits.insert(epochAsCostLimits);
} else {
@@ -398,383 +424,7 @@ namespace storm {
}
return true;
}
-
-
-
-
-
-
-
-
- ///////
-
-
- template<typename ValueType>
- void filterDominatedPoints(std::vector<std::vector<ValueType>>& points, std::vector<storm::solver::OptimizationDirection> const& dirs) {
- std::vector<std::vector<ValueType>> result;
- // Note: this is slow and not inplace but also most likely not performance critical
- for (auto const& p1 : points) {
- bool p1Dominated = false;
- for (auto const& p2 : points) {
- assert(p1.size() == p2.size());
- bool p2DominatesP1 = false;
- for (uint64_t i = 0; i < dirs.size(); ++i) {
- if (storm::solver::minimize(dirs[i]) ? p2[i] <= p1[i] : p2[i] >= p1[i]) {
- if (p2[i] != p1[i]) {
- p2DominatesP1 = true;
- }
- } else {
- p2DominatesP1 = false;
- break;
- }
- }
- if (p2DominatesP1) {
- p1Dominated = true;
- break;
- }
- }
- if (!p1Dominated) {
- result.push_back(p1);
- }
- }
- points = std::move(result);
- }
-
- template<typename ModelType>
- std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment& env) const {
- std::vector<std::vector<ValueType>> result;
-
- auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
- auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
- std::vector<storm::storage::BitVector> dimensionsAsBitVector;
- std::vector<uint64_t> dimensions;
- std::vector<storm::solver::OptimizationDirection> optimizationDirections;
- std::vector<bool> lowerCostBounds;
- for (auto const& dirVar : quantileFormula.getBoundVariables()) {
- dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
- STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
- dimensions.push_back(*dimensionsAsBitVector.back().begin());
- lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
- optimizationDirections.push_back(dirVar.first);
- }
- STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
- if (optimizationDirections[0] == optimizationDirections[1]) {
- if (lowerCostBounds[0] == lowerCostBounds[1]) {
- // TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0])
-
- bool infInfProbSatisfiesFormula = checkLimitValue(env, storm::storage::BitVector(getDimension(), false));
- bool zeroZeroProbSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1]);
- if (infInfProbSatisfiesFormula != zeroZeroProbSatisfiesFormula) {
- std::vector<std::pair<int64_t, typename ModelType::ValueType>> singleQuantileValues;
- for (uint64_t i = 0; i < 2; ++i) {
- // find out whether the bounds B_i = 0 and B_1-i = infinity satisfy the formula
- bool zeroInfSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[1-i]);
- if (zeroInfSatisfiesFormula == storm::solver::minimize(optimizationDirections[0])) {
- // There is bound value b such that the point B_i=0 and B_1-i = b is part of the result
- singleQuantileValues.emplace_back(0, storm::utility::zero<ValueType>());
- } else {
- // Compute quantile where 1-i is set to infinity
- singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]).get());
- if (!storm::solver::minimize(optimizationDirections[i])) {
- // When maximizing bounds, the computed single dimensional quantile value is sat for all values of the other bound.
- // Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i]
- ++singleQuantileValues.back().first;
- }
- }
- // Decrease the value for lower cost bounds to convert >= to >
- if (lowerCostBounds[i]) {
- --singleQuantileValues[i].first;
- }
- }
- std::vector<int64_t> currentEpochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0).
- while (!exploreTwoDimensionalQuantile(env, singleQuantileValues, currentEpochValues, result)) {
- increasePrecision(env);
- ++numPrecisionRefinements;
- }
- } else if (infInfProbSatisfiesFormula) {
- // then also zeroZeroProb satisfies the formula, i.e., all bound values satisfy the formula
- if (storm::solver::minimize(optimizationDirections[0])) {
- result = {{storm::utility::zero<ValueType>(), storm::utility::zero<ValueType>()}};
- } else {
- result = {{ storm::utility::infinity<ValueType>(), storm::utility::infinity<ValueType>()}};
- }
- } else {
- // i.e., no bound value satisfies the formula
- result = {{}};
- }
- } else {
- // TODO: this is an "unreasonable" case
- }
- } else {
- // TODO: find reasonable min/max cases
- }
- return result;
- }
-
- template<typename ModelType>
- bool QuantileHelper<ModelType>::exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const {
- // Initialize some data for easy access
- auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
- auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
- std::vector<storm::storage::BitVector> dimensionsAsBitVector;
- std::vector<uint64_t> dimensions;
- std::vector<storm::solver::OptimizationDirection> optimizationDirections;
- std::vector<bool> lowerCostBounds;
- for (auto const& dirVar : quantileFormula.getBoundVariables()) {
- dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
- STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
- dimensions.push_back(*dimensionsAsBitVector.back().begin());
- lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
- optimizationDirections.push_back(dirVar.first);
- }
- STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
-
-
- MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None)));
-
- // initialize data that will be needed for each epoch
- auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
- auto upperBound = rewardUnfolding.getUpperObjectiveBound();
- std::vector<ValueType> x, b;
- std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
-
- if (currentEpochValues[0] < 0 && currentEpochValues[1] < 0) {
- // This case can only happen in these cases:
- assert(lowerCostBounds[0]);
- assert(currentEpochValues[0] == -1);
- assert(currentEpochValues[1] == -1);
- // This case has been checked already, so we can skip it.
- // Skipping this is actually necessary, since the rewardUnfolding will handle formulas like [F{"a"}>A,{"b"}>B "init"] incorrectly if A=B=-1.
- ++currentEpochValues[1];
- }
- std::set<typename EpochManager::Epoch> exploredEpochs;
- while (true) {
- auto currEpoch = rewardUnfolding.getStartEpoch(true);
- for (uint64_t i = 0; i < 2; ++i) {
- if (currentEpochValues[i] >= 0) {
- rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) currentEpochValues[i]);
- } else {
- rewardUnfolding.getEpochManager().setBottomDimension(currEpoch, dimensions[i]);
- }
- }
- auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
- for (auto const& epoch : epochOrder) {
- if (exploredEpochs.count(epoch) == 0) {
- auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
- ++numCheckedEpochs;
- exploredEpochs.insert(epoch);
- }
- }
- ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
- bool propertySatisfied;
- if (env.solver().isForceSoundness()) {
- auto lowerUpperValue = getLowerUpperBound(env, currValue);
- propertySatisfied = probOpFormula.getBound().isSatisfied(lowerUpperValue.first);
- if (propertySatisfied != probOpFormula.getBound().isSatisfied(lowerUpperValue.second)) {
- // unclear result due to insufficient precision.
- return false;
- }
- } else {
- propertySatisfied = probOpFormula.getBound().isSatisfied(currValue);
- }
-
- // If the reward bounds should be as small as possible, we should stop as soon as the property is satisfied.
- // If the reward bounds should be as large as possible, we should stop as soon as the property is violated and then go a step backwards
- if (storm::solver::minimize(optimizationDirections[0]) == propertySatisfied) {
- // We found another point for the result! Translate it to the original domain
- auto point = currentEpochValues;
- std::vector<ValueType> convertedPoint;
- for (uint64_t i = 0; i < 2; ++i) {
- if (lowerCostBounds[i]) {
- // Translate > to >=
- ++point[i];
- }
- if (i == 1 && storm::solver::maximize(optimizationDirections[i])) {
- // When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction
- --point[i];
- }
- if (point[i] < 0) {
- // Skip this point
- convertedPoint.clear();
- continue;
- }
- convertedPoint.push_back(storm::utility::convertNumber<ValueType>(point[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor);
- }
- if (!convertedPoint.empty()) {
- resultPoints.push_back(std::move(convertedPoint));
- }
-
- if (currentEpochValues[1] == startEpochValues[1].first) {
- break;
- } else {
- ++currentEpochValues[0];
- currentEpochValues[1] = startEpochValues[1].first;
- }
- } else {
- ++currentEpochValues[1];
- }
- }
-
-
- // When maximizing, there are border cases where one dimension can be arbitrarily large
- for (uint64_t i = 0; i < 2; ++i) {
- if (storm::solver::maximize(optimizationDirections[i]) && startEpochValues[i].first > 0) {
- std::vector<ValueType> newResultPoint(2);
- newResultPoint[i] = storm::utility::convertNumber<ValueType>(startEpochValues[i].first - 1) * startEpochValues[i].second;
- newResultPoint[1-i] = storm::utility::infinity<ValueType>();
- resultPoints.push_back(newResultPoint);
- }
- }
-
- filterDominatedPoints(resultPoints, optimizationDirections);
- return true;
- }
- template<typename ModelType>
- boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantileForDimension(Environment const& env, uint64_t dimension) const {
- // We assume that there is one bound value that violates the quantile and one bound value that satisfies it.
-
- // Let all other open bounds approach infinity
- std::set<storm::expressions::Variable> infinityVariables;
- for (auto const& d : getOpenDimensions()) {
- if (d != dimension) {
- infinityVariables.insert(getVariableForDimension(d));
- }
- }
- auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None));
- MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
-
- bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
- bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
-
- // initialize data that will be needed for each epoch
- auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
- auto upperBound = rewardUnfolding.getUpperObjectiveBound();
- std::vector<ValueType> x, b;
- std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
-
- uint64_t epochValue = 0;
- std::set<typename EpochManager::Epoch> exploredEpochs;
- while (true) {
- auto currEpoch = rewardUnfolding.getStartEpoch(true);
- rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimension, (uint64_t) epochValue);
- auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
- for (auto const& epoch : epochOrder) {
- if (exploredEpochs.count(epoch) == 0) {
- auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
- ++numCheckedEpochs;
- exploredEpochs.insert(epoch);
- }
- }
- ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
-
- bool propertySatisfied;
- if (env.solver().isForceSoundness()) {
- auto lowerUpperValue = getLowerUpperBound(env, currValue);
- propertySatisfied = transformedFormula->getBound().isSatisfied(lowerUpperValue.first);
- if (propertySatisfied != transformedFormula->getBound().isSatisfied(lowerUpperValue.second)) {
- // unclear result due to insufficient precision.
- return boost::none;
- }
- } else {
- propertySatisfied = transformedFormula->getBound().isSatisfied(currValue);
- }
- // If the reward bound should be as small as possible, we should stop as soon as the property is satisfied.
- // If the reward bound should be as large as possible, we should stop as soon as sthe property is violated and then go a step backwards
- if (minimizingRewardBound && propertySatisfied) {
- // We found a satisfying value!
- if (!upperCostBound) {
- // The rewardunfolding assumes strict lower cost bounds while we assume non-strict ones. Hence, >B becomes >=(B+1).
- ++epochValue;
- }
- break;
- } else if (!minimizingRewardBound && !propertySatisfied) {
- // We found a non-satisfying value. Go one step back to get the largest satisfying value.
- // ... however, lower cost bounds need to be converted from strict bounds >B to non strict bounds >=(B+1).
- // Hence, no operation is necessary in this case.
- if (upperCostBound) {
- STORM_LOG_ASSERT(epochValue > 0, "The property does not seem to be satisfiable. This case should have been excluded earlier.");
- --epochValue;
- }
- break;
- }
- ++epochValue;
- }
- return std::make_pair(epochValue, rewardUnfolding.getDimension(dimension).scalingFactor);
- }
-
- template<typename ModelType>
- bool QuantileHelper<ModelType>::checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const {
- auto const& probabilityBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound();
- // Increase the precision until we get a conclusive result
- while (true) {
- ValueType numericResult = computeLimitValue(env, infDimensions);
- if (env.solver().isForceSoundness()) {
- auto lowerUpper = getLowerUpperBound(env, numericResult);
- bool lowerSat = probabilityBound.isSatisfied(lowerUpper.first);
- bool upperSat = probabilityBound.isSatisfied(lowerUpper.second);
- if (lowerSat == upperSat) {
- return lowerSat;
- } else {
- increasePrecision(env);
- ++numPrecisionRefinements;
- }
- } else {
- return probabilityBound.isSatisfied(numericResult);
- }
- }
- }
-
- template<typename ModelType>
- typename ModelType::ValueType QuantileHelper<ModelType>::computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const {
- // To compute the limit for an upper bounded dimension, we can simply drop the bound
- // To compute the limit for a lower bounded dimension, we only consider epoch steps that lie in an end component and check for the bound >0 instead.
- // Notice, however, that this approach becomes problematic if, at the same time, we consider an upper bounded dimension with bound value zero.
- std::vector<BoundTransformation> bts(getDimension(), BoundTransformation::None);
- std::set<storm::expressions::Variable> infinityVariables;
- for (auto const& d : getOpenDimensions()) {
- bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d);
- if (infDimensions.get(d)) {
- infinityVariables.insert(getVariableForDimension(d));
- } else {
- if (upperCostBound) {
- bts[d] = BoundTransformation::LessEqualZero;
- } else {
- bts[d] = BoundTransformation::GreaterEqualZero;
- }
- }
- }
- auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts);
-
- MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
- if (!rewardUnfolding.getProb1Objectives().empty()) {
- assert(rewardUnfolding.getProb1Objectives().size() == 1);
- // The probability is one.
- return storm::utility::one<ValueType>();
- }
- // Get lower and upper bounds for the solution.
- auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
- auto upperBound = rewardUnfolding.getUpperObjectiveBound();
-
- // Initialize epoch models
- auto initEpoch = rewardUnfolding.getStartEpoch();
- auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
- // initialize data that will be needed for each epoch
- std::vector<ValueType> x, b;
- std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
-
- for (auto const& epoch : epochOrder) {
- auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
- ++numCheckedEpochs;
- }
-
- ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch);
- STORM_LOG_TRACE("Limit probability for infinity dimensions " << infDimensions << " is " << numericResult << ".");
- return numericResult;
- }
template class QuantileHelper<storm::models::sparse::Mdp<double>>;
template class QuantileHelper<storm::models::sparse::Mdp<storm::RationalNumber>>;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index 7b8b254cc..f9f92163c 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -1,17 +1,17 @@
#pragma once
+#include <map>
+#include <boost/optional.hpp>
+
#include "storm/logic/QuantileFormula.h"
-#include "boost/optional.hpp"
+#include "storm/logic/ProbabilityOperatorFormula.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
+#include "storm/storage/BitVector.h"
namespace storm {
class Environment;
- namespace storage {
- class BitVector;
- }
-
namespace modelchecker {
namespace helper {
namespace rewardbounded {
@@ -22,12 +22,12 @@ namespace storm {
public:
QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula);
- std::vector<std::vector<ValueType>> computeQuantile(Environment const& env) const;
+ std::vector<std::vector<ValueType>> computeQuantile(Environment const& env);
private:
- std::pair<CostLimitClosure, std::vector<ValueType>> computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions) const;
- bool computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding);
+ std::pair<CostLimitClosure, std::vector<ValueType>> computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, bool complementaryQuery);
+ bool computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding);
std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env) const;
@@ -66,6 +66,7 @@ namespace storm {
ModelType const& model;
storm::logic::QuantileFormula const& quantileFormula;
+ std::map<storm::storage::BitVector, std::pair<CostLimitClosure, std::vector<ValueType>>> cachedSubQueryResults;
/// Statistics
mutable uint64_t numCheckedEpochs;
From c33ac18a5a63e391120bb1365b58130305e4916f Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Wed, 13 Mar 2019 09:35:07 +0100
Subject: [PATCH 03/11] Quantiles: Fixed a precision related issue in new
implementation.
---
.../prctl/helper/rewardbounded/EpochManager.cpp | 11 +++++++++++
.../prctl/helper/rewardbounded/EpochManager.h | 1 +
.../MultiDimensionalRewardUnfolding.cpp | 6 +-----
.../helper/rewardbounded/QuantileHelper.cpp | 16 ++++++++++------
4 files changed, 23 insertions(+), 11 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
index e248a74da..3537b618e 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
@@ -210,6 +210,17 @@ namespace storm {
return (epoch >> (dimension * bitsPerDimension)) & dimensionBitMask;
}
+ uint64_t EpochManager::getSumOfDimensions(Epoch const& epoch) const {
+ STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+ uint64_t sumOfDimensions = 0;
+ for (uint64_t dim = 0; dim < getDimensionCount(); ++dim) {
+ if (!isBottomDimension(epoch, dim)) {
+ sumOfDimensions += getDimensionOfEpoch(epoch, dim) + 1;
+ }
+ }
+ return sumOfDimensions;
+ }
+
std::string EpochManager::toString(Epoch const& epoch) const {
STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
std::string res = "<" + (isBottomDimension(epoch, 0) ? "_" : std::to_string(getDimensionOfEpoch(epoch, 0)));
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
index 6784cdac4..b1ff3e5ce 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
@@ -45,6 +45,7 @@ namespace storm {
bool isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const;
bool isBottomDimensionEpochClass(EpochClass const& epochClass, uint64_t const& dimension) const;
uint64_t getDimensionOfEpoch(Epoch const& epoch, uint64_t const& dimension) const; // assumes that the dimension is not bottom
+ uint64_t getSumOfDimensions(Epoch const& epoch) const; // assumes that the dimension is not bottom
std::string toString(Epoch const& epoch) const;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
index 5d82672eb..fad4feda2 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
@@ -711,11 +711,7 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
ValueType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) {
- uint64_t sumOfDimensions = 0;
- for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
- sumOfDimensions += epochManager.getDimensionOfEpoch(startEpoch, dim) + 1;
- }
- return precision / storm::utility::convertNumber<ValueType>(sumOfDimensions);
+ return precision / storm::utility::convertNumber<ValueType>(epochManager.getSumOfDimensions(startEpoch) + 1);
}
template<typename ValueType, bool SingleObjectiveMode>
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index e78354752..d2f6c1c80 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -64,9 +64,9 @@ namespace storm {
}
// TODO
- // Fix precision increasing
// Multiple quantile formulas in the same file yield constants def clash
// ignore optimization direction for quantile variables
+ // Test cases
}
enum class BoundTransformation {
@@ -139,9 +139,12 @@ namespace storm {
env.solver().minMax().setPrecision(env.solver().minMax().getPrecision() * factor);
}
- /// Computes a lower/ upper bound on the actual result of a minmax or linear equation solver
+ /*!
+ * Computes a lower / upper bound on the actual result of a sound minmax or linear equation solver
+ *
+ */
template<typename ValueType>
- std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& value, bool minMax = true) {
+ std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& factor, ValueType const& value, bool minMax = true) {
ValueType prec;
bool relative;
if (minMax) {
@@ -152,9 +155,9 @@ namespace storm {
relative = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).second.get();
}
if (relative) {
- return std::make_pair<ValueType, ValueType>(value * (1/(prec + 1)), value * (1 + prec/(prec +1)));
+ return std::make_pair<ValueType, ValueType>((value * (1/(prec + 1))) * factor, (value * (1 + prec/(prec +1))) * factor);
} else {
- return std::make_pair<ValueType, ValueType>(value - prec, value + prec);
+ return std::make_pair<ValueType, ValueType>((value - prec) * factor, (value + prec) * factor);
}
}
@@ -402,7 +405,8 @@ namespace storm {
ValueType currValue = rewardUnfolding.getInitialStateResult(epoch);
bool propertySatisfied;
if (env.solver().isForceSoundness()) {
- auto lowerUpperValue = getLowerUpperBound(env, currValue);
+ ValueType sumOfEpochDimensions = storm::utility::convertNumber<ValueType>(rewardUnfolding.getEpochManager().getSumOfDimensions(epoch) + 1);
+ auto lowerUpperValue = getLowerUpperBound(env, sumOfEpochDimensions, currValue);
propertySatisfied = boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.first);
if (propertySatisfied != boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.second)) {
// unclear result due to insufficient precision.
From 38121c28cb0a2acb053216827a10e37c28fad693 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Wed, 13 Mar 2019 18:02:13 +0100
Subject: [PATCH 04/11] quantiles: permute point entries if the order of
quantile variable definitions is not the same as the order of occurrence on a
cost bound.
---
.../helper/rewardbounded/QuantileHelper.cpp | 24 ++++++++++++-------
1 file changed, 16 insertions(+), 8 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index d2f6c1c80..7fc19c5ab 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -208,15 +208,23 @@ namespace storm {
// Call the internal recursive function
auto internalResult = computeQuantile(envCpy, getOpenDimensions(), false);
- // Translate the result by applying the scaling factors.
+ // Translate the result by applying the scaling factors and permutation.
+ std::vector<uint64_t> permutation;
+ for (auto const& v : quantileFormula.getBoundVariables()) {
+ for (auto const& dim : getOpenDimensions()) {
+ if (getVariableForDimension(dim) == v) {
+ permutation.push_back(dim);
+ break;
+ }
+ }
+ }
+ assert(permutation.size() == getOpenDimensions().getNumberOfSetBits());
for (auto const& costLimits : internalResult.first.getGenerator()) {
- std::vector<ValueType> resultPoint(costLimits.size());
- storm::utility::vector::applyPointwise<CostLimit, ValueType, ValueType>(costLimits, internalResult.second, resultPoint, [](CostLimit const& costLimit, ValueType const& factor) -> ValueType {
- if (costLimit.isInfinity()) {
- return storm::utility::infinity<ValueType>();
- } else {
- return storm::utility::convertNumber<ValueType>(costLimit.get()) * factor;
- }});
+ std::vector<ValueType> resultPoint;
+ for (auto const& dim : permutation) {
+ CostLimit const& cl = costLimits[dim];
+ resultPoint.push_back(cl.isInfinity() ? storm::utility::infinity<ValueType>() : storm::utility::convertNumber<ValueType>(cl.get()) * internalResult.second[dim]);
+ }
result.push_back(resultPoint);
}
if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
From 8a72aee7648d06690e5a6f916dea5e0c47932866 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Wed, 13 Mar 2019 18:04:16 +0100
Subject: [PATCH 05/11] QuantileFormulas: ignore optimization direction
(min/max) for quantile variables.
---
.../parser/FormulaParserGrammar.cpp | 15 +++++++-------
.../parser/FormulaParserGrammar.h | 6 +++---
src/storm/logic/QuantileFormula.cpp | 20 +++++--------------
src/storm/logic/QuantileFormula.h | 9 +++------
.../helper/rewardbounded/QuantileHelper.cpp | 20 +++----------------
.../helper/rewardbounded/QuantileHelper.h | 20 -------------------
6 files changed, 21 insertions(+), 69 deletions(-)
diff --git a/src/storm-parsers/parser/FormulaParserGrammar.cpp b/src/storm-parsers/parser/FormulaParserGrammar.cpp
index 6f565fc20..d2db73a6a 100644
--- a/src/storm-parsers/parser/FormulaParserGrammar.cpp
+++ b/src/storm-parsers/parser/FormulaParserGrammar.cpp
@@ -127,9 +127,9 @@ namespace storm {
identifier %= qi::as_string[qi::raw[qi::lexeme[((qi::alpha | qi::char_('_') | qi::char_('.')) >> *(qi::alnum | qi::char_('_')))]]];
identifier.name("identifier");
- quantileBoundVariable = ((qi::lit("min")[qi::_a = storm::solver::OptimizationDirection::Minimize] | qi::lit("max")[qi::_a = storm::solver::OptimizationDirection::Maximize]) >> identifier)[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileBoundVariables, phoenix::ref(*this), qi::_a, qi::_1)];
- quantileBoundVariable.name("Quantile bound variable");
- quantileFormula = (qi::lit("quantile") > qi::lit("(") >> (quantileBoundVariable % qi::lit(",")) >> qi::lit(",") >> stateFormula > qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileFormula, phoenix::ref(*this), qi::_1, qi::_2)];
+ quantileBoundVariable = (-(qi::lit("min")[qi::_a = storm::solver::OptimizationDirection::Minimize] | qi::lit("max")[qi::_a = storm::solver::OptimizationDirection::Maximize]) >> identifier >> qi::lit(","))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileBoundVariables, phoenix::ref(*this), qi::_a, qi::_1)];
+ quantileBoundVariable.name("quantile bound variable");
+ quantileFormula = (qi::lit("quantile") > qi::lit("(") >> *(quantileBoundVariable) >> stateFormula > qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileFormula, phoenix::ref(*this), qi::_1, qi::_2)];
quantileFormula.name("Quantile formula");
stateFormula = (orStateFormula | multiFormula | quantileFormula);
@@ -424,17 +424,16 @@ namespace storm {
}
}
- std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> FormulaParserGrammar::createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName) {
+ storm::expressions::Variable FormulaParserGrammar::createQuantileBoundVariables(boost::optional<storm::solver::OptimizationDirection> const& dir, std::string const& variableName) {
STORM_LOG_ASSERT(manager, "Mutable expression manager required to define quantile bound variable.");
STORM_LOG_THROW(!manager->hasVariable(variableName), storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in property: variable already exists.");
-
+ STORM_LOG_WARN_COND(!dir.is_initialized(), "Optimization direction '" << dir.get() << "' for quantile variable " << variableName << " is ignored. This information will be derived from the subformula of the quantile.");
storm::expressions::Variable var = manager->declareRationalVariable(variableName);
addIdentifierExpression(variableName, var);
- std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> result(dir, var);
- return result;
+ return var;
}
- std::shared_ptr<storm::logic::Formula const> FormulaParserGrammar::createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula) {
+ std::shared_ptr<storm::logic::Formula const> FormulaParserGrammar::createQuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula) {
return std::shared_ptr<storm::logic::Formula const>(new storm::logic::QuantileFormula(boundVariables, subformula));
}
diff --git a/src/storm-parsers/parser/FormulaParserGrammar.h b/src/storm-parsers/parser/FormulaParserGrammar.h
index e582abe40..7f2bbf33c 100644
--- a/src/storm-parsers/parser/FormulaParserGrammar.h
+++ b/src/storm-parsers/parser/FormulaParserGrammar.h
@@ -197,7 +197,7 @@ namespace storm {
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> longRunAverageRewardFormula;
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> multiFormula;
- qi::rule<Iterator, std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>(), qi::locals<storm::solver::OptimizationDirection>, Skipper> quantileBoundVariable;
+ qi::rule<Iterator, storm::expressions::Variable(), qi::locals<boost::optional<storm::solver::OptimizationDirection>>, Skipper> quantileBoundVariable;
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> quantileFormula;
// Parser that is used to recognize doubles only (as opposed to Spirit's double_ parser).
@@ -232,8 +232,8 @@ namespace storm {
std::shared_ptr<storm::logic::Formula const> createBinaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& leftSubformula, std::shared_ptr<storm::logic::Formula const> const& rightSubformula, storm::logic::BinaryBooleanStateFormula::OperatorType operatorType);
std::shared_ptr<storm::logic::Formula const> createUnaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& subformula, boost::optional<storm::logic::UnaryBooleanStateFormula::OperatorType> const& operatorType);
std::shared_ptr<storm::logic::Formula const> createMultiFormula(std::vector<std::shared_ptr<storm::logic::Formula const>> const& subformulas);
- std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName);
- std::shared_ptr<storm::logic::Formula const> createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula);
+ storm::expressions::Variable createQuantileBoundVariables(boost::optional<storm::solver::OptimizationDirection> const& dir, std::string const& variableName);
+ std::shared_ptr<storm::logic::Formula const> createQuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula);
std::set<storm::expressions::Variable> getUndefinedConstants(std::shared_ptr<storm::logic::Formula const> const& formula) const;
storm::jani::Property createProperty(boost::optional<std::string> const& propertyName, storm::modelchecker::FilterType const& filterType, std::shared_ptr<storm::logic::Formula const> const& formula, std::shared_ptr<storm::logic::Formula const> const& states);
diff --git a/src/storm/logic/QuantileFormula.cpp b/src/storm/logic/QuantileFormula.cpp
index 322f36fdd..4e1e21aef 100644
--- a/src/storm/logic/QuantileFormula.cpp
+++ b/src/storm/logic/QuantileFormula.cpp
@@ -7,7 +7,7 @@
namespace storm {
namespace logic {
- QuantileFormula::QuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> subformula) : boundVariables(boundVariables), subformula(subformula) {
+ QuantileFormula::QuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<Formula const> subformula) : boundVariables(boundVariables), subformula(subformula) {
STORM_LOG_THROW(!boundVariables.empty(), storm::exceptions::InvalidArgumentException, "Quantile formula without bound variables are invalid.");
}
@@ -45,28 +45,18 @@ namespace storm {
storm::expressions::Variable const& QuantileFormula::getBoundVariable() const {
STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique bound variables. However, there are multiple bound variables defined.");
- return boundVariables.front().second;
+ return boundVariables.front();
}
storm::expressions::Variable const& QuantileFormula::getBoundVariable(uint64_t index) const {
STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested bound variable with index" << index << ". However, there are only " << boundVariables.size() << " bound variables.");
- return boundVariables[index].second;
+ return boundVariables[index];
}
- std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& QuantileFormula::getBoundVariables() const {
+ std::vector<storm::expressions::Variable> const& QuantileFormula::getBoundVariables() const {
return boundVariables;
}
- storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection() const {
- STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique optimization direction of the bound variables. However, there are multiple bound variables defined.");
- return boundVariables.front().first;
- }
-
- storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection(uint64_t index) const {
- STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested optimization direction with index" << index << ". However, there are only " << boundVariables.size() << " bound variables.");
- return boundVariables[index].first;
- }
-
boost::any QuantileFormula::accept(FormulaVisitor const& visitor, boost::any const& data) const {
return visitor.visit(*this, data);
}
@@ -86,7 +76,7 @@ namespace storm {
std::ostream& QuantileFormula::writeToStream(std::ostream& out) const {
out << "quantile(";
for (auto const& bv : boundVariables) {
- out << (storm::solver::minimize(bv.first) ? "min" : "max") << " " << bv.second.getName() << ", ";
+ out << bv.getName() << ", ";
}
subformula->writeToStream(out);
out << ")";
diff --git a/src/storm/logic/QuantileFormula.h b/src/storm/logic/QuantileFormula.h
index 2690171a5..8765013a7 100644
--- a/src/storm/logic/QuantileFormula.h
+++ b/src/storm/logic/QuantileFormula.h
@@ -1,13 +1,12 @@
#pragma once
#include "storm/logic/Formula.h"
-#include "storm/solver/OptimizationDirection.h"
namespace storm {
namespace logic {
class QuantileFormula : public Formula {
public:
- QuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> subformula);
+ QuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<Formula const> subformula);
virtual ~QuantileFormula();
@@ -23,9 +22,7 @@ namespace storm {
storm::expressions::Variable const& getBoundVariable() const;
storm::expressions::Variable const& getBoundVariable(uint64_t index) const;
- std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& getBoundVariables() const;
- storm::solver::OptimizationDirection const& getOptimizationDirection() const;
- storm::solver::OptimizationDirection const& getOptimizationDirection(uint64_t index) const;
+ std::vector<storm::expressions::Variable> const& getBoundVariables() const;
virtual boost::any accept(FormulaVisitor const& visitor, boost::any const& data) const override;
virtual void gatherAtomicExpressionFormulas(std::vector<std::shared_ptr<AtomicExpressionFormula const>>& atomicExpressionFormulas) const override;
@@ -34,7 +31,7 @@ namespace storm {
virtual std::ostream& writeToStream(std::ostream& out) const override;
private:
- std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> boundVariables;
+ std::vector<storm::expressions::Variable> boundVariables;
std::shared_ptr<Formula const> subformula;
};
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index 7fc19c5ab..f965d9794 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -32,15 +32,14 @@ namespace storm {
// Do all kinds of sanity check.
std::set<storm::expressions::Variable> quantileVariables;
for (auto const& quantileVariable : quantileFormula.getBoundVariables()) {
- STORM_LOG_THROW(quantileVariables.count(quantileVariable.second) == 0, storm::exceptions::NotSupportedException, "Quantile formula considers the same bound variable twice.");
- quantileVariables.insert(quantileVariable.second);
+ STORM_LOG_THROW(quantileVariables.count(quantileVariable) == 0, storm::exceptions::NotSupportedException, "Quantile formula considers the same bound variable twice.");
+ quantileVariables.insert(quantileVariable);
}
STORM_LOG_THROW(quantileFormula.getSubformula().isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs probability operator inside. The formula " << quantileFormula << " is not supported.");
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
STORM_LOG_THROW(probOpFormula.hasBound(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have a bound.");
STORM_LOG_THROW(!model.isNondeterministicModel() || probOpFormula.hasOptimalityType(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have an optimality type.");
STORM_LOG_WARN_COND(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, "Probability operator inside quantile formula needs to have bound > or <=. The specified comparison type might lead to non-termination."); // This has to do with letting bound variables approach infinity, e.g., Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
- bool lowerBounded = storm::logic::isLowerBound(probOpFormula.getBound().comparisonType);
STORM_LOG_THROW(probOpFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs bounded until probability operator formula as subformula. The formula " << quantileFormula << " is not supported.");
auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
std::set<storm::expressions::Variable> boundVariables;
@@ -65,7 +64,6 @@ namespace storm {
// TODO
// Multiple quantile formulas in the same file yield constants def clash
- // ignore optimization direction for quantile variables
// Test cases
}
@@ -178,19 +176,7 @@ namespace storm {
}
return res;
}
-
- template<typename ModelType>
- storm::solver::OptimizationDirection const& QuantileHelper<ModelType>::getOptimizationDirForDimension(uint64_t const& dim) const {
- storm::expressions::Variable const& dimVar = getVariableForDimension(dim);
- for (auto const& boundVar : quantileFormula.getBoundVariables()) {
- if (boundVar.second == dimVar) {
- return boundVar.first;
- }
- }
- STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "The bound variable '" << dimVar.getName() << "' is not specified within the quantile formula '" << quantileFormula << "'.");
- return quantileFormula.getOptimizationDirection();
- }
-
+
template<typename ModelType>
storm::expressions::Variable const& QuantileHelper<ModelType>::getVariableForDimension(uint64_t const& dim) const {
auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index f9f92163c..e3bb5fd78 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -30,25 +30,6 @@ namespace storm {
bool computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding);
- std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env) const;
- bool exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const;
-
- /*!
- * Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
- */
- ValueType computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const;
-
- /*!
- * Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
- * The computed value is compared to the probability threshold.
- * In sound mode, precision is iteratively increased in case of 'inconsistent' results.
- */
- bool checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const;
-
- /// Computes the quantile with respect to the given dimension.
- /// boost::none is returned in case of insufficient precision.
- boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> computeQuantileForDimension(Environment const& env, uint64_t dim) const;
-
/*!
* Gets the number of dimensions of the underlying boudned until formula
*/
@@ -61,7 +42,6 @@ namespace storm {
storm::storage::BitVector getOpenDimensions() const;
storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var) const;
- storm::solver::OptimizationDirection const& getOptimizationDirForDimension(uint64_t const& dim) const;
storm::expressions::Variable const& getVariableForDimension(uint64_t const& dim) const;
ModelType const& model;
From c21ea2ce1f54a4811169325b0220559e3e0c1b16 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Wed, 13 Mar 2019 18:36:26 +0100
Subject: [PATCH 06/11] Quantiles: Bug fixes.
---
.../prctl/helper/rewardbounded/QuantileHelper.cpp | 9 ++++++---
1 file changed, 6 insertions(+), 3 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index f965d9794..eefd43d95 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -152,10 +152,11 @@ namespace storm {
prec = storm::utility::convertNumber<ValueType>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get());
relative = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).second.get();
}
+ prec *= factor;
if (relative) {
- return std::make_pair<ValueType, ValueType>((value * (1/(prec + 1))) * factor, (value * (1 + prec/(prec +1))) * factor);
+ return std::make_pair<ValueType, ValueType>(value * (1/(prec + 1)), value * (1 + prec/(prec +1)));
} else {
- return std::make_pair<ValueType, ValueType>((value - prec) * factor, (value + prec) * factor);
+ return std::make_pair<ValueType, ValueType>(value - prec, value + prec);
}
}
@@ -197,11 +198,13 @@ namespace storm {
// Translate the result by applying the scaling factors and permutation.
std::vector<uint64_t> permutation;
for (auto const& v : quantileFormula.getBoundVariables()) {
+ uint64_t openDim = 0;
for (auto const& dim : getOpenDimensions()) {
if (getVariableForDimension(dim) == v) {
- permutation.push_back(dim);
+ permutation.push_back(openDim);
break;
}
+ ++openDim;
}
}
assert(permutation.size() == getOpenDimensions().getNumberOfSetBits());
From 971f4c85084717b2a1d556f0105109aaeaef6b67 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Thu, 14 Mar 2019 12:36:20 +0100
Subject: [PATCH 07/11] Quantiles: Fixed analysing epochs unnecessarily, fixed
having multiple quantile formulas over the same variables.
---
.../parser/FormulaParserGrammar.cpp | 10 +++-
.../parser/FormulaParserGrammar.h | 1 +
.../helper/rewardbounded/CostLimitClosure.cpp | 14 ++++-
.../helper/rewardbounded/CostLimitClosure.h | 2 +
.../helper/rewardbounded/QuantileHelper.cpp | 60 ++++++++++---------
.../helper/rewardbounded/QuantileHelper.h | 2 +
6 files changed, 57 insertions(+), 32 deletions(-)
diff --git a/src/storm-parsers/parser/FormulaParserGrammar.cpp b/src/storm-parsers/parser/FormulaParserGrammar.cpp
index d2db73a6a..69c317a21 100644
--- a/src/storm-parsers/parser/FormulaParserGrammar.cpp
+++ b/src/storm-parsers/parser/FormulaParserGrammar.cpp
@@ -426,9 +426,15 @@ namespace storm {
storm::expressions::Variable FormulaParserGrammar::createQuantileBoundVariables(boost::optional<storm::solver::OptimizationDirection> const& dir, std::string const& variableName) {
STORM_LOG_ASSERT(manager, "Mutable expression manager required to define quantile bound variable.");
- STORM_LOG_THROW(!manager->hasVariable(variableName), storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in property: variable already exists.");
+ storm::expressions::Variable var;
+ if (manager->hasVariable(variableName)) {
+ var = manager->getVariable(variableName);
+ STORM_LOG_THROW(quantileFormulaVariables.count(var) > 0, storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in quantile formula: variable already exists.");
+ } else {
+ var = manager->declareRationalVariable(variableName);
+ quantileFormulaVariables.insert(var);
+ }
STORM_LOG_WARN_COND(!dir.is_initialized(), "Optimization direction '" << dir.get() << "' for quantile variable " << variableName << " is ignored. This information will be derived from the subformula of the quantile.");
- storm::expressions::Variable var = manager->declareRationalVariable(variableName);
addIdentifierExpression(variableName, var);
return var;
}
diff --git a/src/storm-parsers/parser/FormulaParserGrammar.h b/src/storm-parsers/parser/FormulaParserGrammar.h
index 7f2bbf33c..64a38e4c1 100644
--- a/src/storm-parsers/parser/FormulaParserGrammar.h
+++ b/src/storm-parsers/parser/FormulaParserGrammar.h
@@ -245,6 +245,7 @@ namespace storm {
uint64_t propertyCount;
std::set<storm::expressions::Variable> undefinedConstants;
+ std::set<storm::expressions::Variable> quantileFormulaVariables;
};
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
index f4fd31d85..fefc20299 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
@@ -95,7 +95,19 @@ namespace storm {
}
return contains(infinityProjection);
}
-
+
+ bool CostLimitClosure::empty() const {
+ return generator.empty();
+ }
+
+ bool CostLimitClosure::full() const {
+ CostLimits p(dimension(), CostLimit(0));
+ for (auto const& dim : downwardDimensions) {
+ p[dim] = CostLimit::infinity();
+ }
+ return contains(p);
+ }
+
bool CostLimitClosure::dominates(CostLimits const& lhs, CostLimits const& rhs) const {
for (uint64_t i = 0; i < lhs.size(); ++i) {
if (downwardDimensions.get(i)) {
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
index 579a33b5b..06abb4722 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
@@ -38,6 +38,8 @@ namespace storm {
bool contains(CostLimits const& costLimits) const;
bool containsUpwardClosure(CostLimits const& costLimits) const;
bool dominates(CostLimits const& lhs, CostLimits const& rhs) const;
+ bool empty() const; // True if there is no point in this
+ bool full() const; // True if all points lie in this.
std::vector<CostLimits> getDominatingCostLimits(CostLimits const& costLimits) const;
GeneratorType const& getGenerator() const;
uint64_t dimension() const;
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index eefd43d95..c6bdda24e 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -21,6 +21,7 @@
#include "storm/logic/BoundedUntilFormula.h"
#include "storm/exceptions/NotSupportedException.h"
+#include "storm/exceptions/UnexpectedException.h"
namespace storm {
namespace modelchecker {
@@ -61,10 +62,6 @@ namespace storm {
STORM_LOG_THROW(quantileVariables.count(boundVariable) == 1, storm::exceptions::NotSupportedException, "The formula contains undefined constant '" << boundExpression << "'.");
}
}
-
- // TODO
- // Multiple quantile formulas in the same file yield constants def clash
- // Test cases
}
enum class BoundTransformation {
@@ -188,6 +185,7 @@ namespace storm {
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) {
numCheckedEpochs = 0;
numPrecisionRefinements = 0;
+ swEpochAnalysis.reset();
cachedSubQueryResults.clear();
std::vector<std::vector<ValueType>> result;
@@ -219,6 +217,7 @@ namespace storm {
if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl;
std::cout << "Number of required precision refinements: " << numPrecisionRefinements << std::endl;
+ std::cout << "Time for epoch model analysis: " << swEpochAnalysis << " seconds." << std::endl;
}
return result;
}
@@ -302,30 +301,26 @@ namespace storm {
}
}
- bool getNextCandidateCostLimit(CostLimit const& maxCostLimit, CostLimits& current) {
- if (maxCostLimit.get() == 0) {
+ bool getNextCandidateCostLimit(CostLimit const& candidateCostLimitSum, CostLimits& current) {
+ if (current.size() == 0) {
return false;
}
- storm::storage::BitVector nonMaxEntries = storm::utility::vector::filter<CostLimit>(current, [&maxCostLimit] (CostLimit const& value) -> bool { return value < maxCostLimit; });
- bool allZero = true;
- for (auto const& entry : nonMaxEntries) {
- if (current[entry].get() > 0) {
- --current[entry].get();
- allZero = false;
- break;
- } else {
- current[entry] = CostLimit(maxCostLimit.get() - 1);
- }
+ uint64_t iSum = current.front().get();
+ if (iSum == candidateCostLimitSum.get()) {
+ return false;
}
- if (allZero) {
- nonMaxEntries.increment();
- if (nonMaxEntries.full()) {
- return false;
+ for (uint64_t i = 1; i < current.size(); ++i) {
+ iSum += current[i].get();
+ if (iSum == candidateCostLimitSum.get()) {
+ ++current[i-1].get();
+ uint64_t newVal = current[i].get() - 1;
+ current[i].get() = 0;
+ current.back().get() = newVal;
+ return true;
}
- current = CostLimits(current.size(), maxCostLimit);
- storm::utility::vector::setVectorValues(current, nonMaxEntries, CostLimit(maxCostLimit.get() - 1));
}
- return true;
+ STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The entries of the current cost limit candidate do not sum up to the current candidate sum.");
+ return false;
}
bool translateEpochToCostLimits(EpochManager::Epoch const& epoch, EpochManager::Epoch const& startEpoch,storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, EpochManager const& epochManager, CostLimits& epochAsCostLimits) {
@@ -367,10 +362,14 @@ namespace storm {
std::set<EpochManager::Epoch> checkedEpochs;
bool progress = true;
- for (CostLimit currentMaxCostLimit(0); progress; ++currentMaxCostLimit.get()) {
- CostLimits currentCandidate(satCostLimits.dimension(), currentMaxCostLimit);
- // We can only stop the exploration, if the upward closure of the point in the 'top right corner' is contained in the (un)satCostlLimits.
- progress = !satCostLimits.containsUpwardClosure(currentCandidate) && !unsatCostLimits.containsUpwardClosure(currentCandidate);
+ for (CostLimit candidateCostLimitSum(0); progress; ++candidateCostLimitSum.get()) {
+ CostLimits currentCandidate(satCostLimits.dimension(), CostLimit(0));
+ if (!currentCandidate.empty()) {
+ currentCandidate.back() = candidateCostLimitSum;
+ }
+ // We can still have progress if one of the closures is empty and the other is not full.
+ // This ensures that we do not terminate too early in case that the (un)satCostLimits are initially non-empty.
+ progress = (satCostLimits.empty() && !unsatCostLimits.full()) || (unsatCostLimits.empty() && !satCostLimits.full());
do {
if (!satCostLimits.contains(currentCandidate) && !unsatCostLimits.contains(currentCandidate)) {
progress = true;
@@ -389,13 +388,16 @@ namespace storm {
}
++costLimitIt;
}
+ STORM_LOG_DEBUG("Checking start epoch " << rewardUnfolding.getEpochManager().toString(startEpoch) << ".");
auto epochSequence = rewardUnfolding.getEpochComputationOrder(startEpoch);
for (auto const& epoch : epochSequence) {
if (checkedEpochs.count(epoch) == 0) {
+ ++numCheckedEpochs;
+ swEpochAnalysis.start();
checkedEpochs.insert(epoch);
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env,boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
- ++numCheckedEpochs;
+ swEpochAnalysis.stop();
CostLimits epochAsCostLimits;
if (translateEpochToCostLimits(epoch, startEpoch, consideredDimensions, lowerBoundedDimensions, rewardUnfolding.getEpochManager(), epochAsCostLimits)) {
@@ -421,7 +423,7 @@ namespace storm {
}
}
}
- } while (getNextCandidateCostLimit(currentMaxCostLimit, currentCandidate));
+ } while (getNextCandidateCostLimit(candidateCostLimitSum, currentCandidate));
}
return true;
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index e3bb5fd78..63a06e34e 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -8,6 +8,7 @@
#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
#include "storm/storage/BitVector.h"
+#include "storm/utility/Stopwatch.h"
namespace storm {
class Environment;
@@ -51,6 +52,7 @@ namespace storm {
/// Statistics
mutable uint64_t numCheckedEpochs;
mutable uint64_t numPrecisionRefinements;
+ mutable storm::utility::Stopwatch swEpochAnalysis;
};
}
}
From aa3a1f5ff7cfbe42bf54df11a48acb612c51a4d5 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Thu, 14 Mar 2019 13:42:23 +0100
Subject: [PATCH 08/11] Quantiles: Improved performance by excluding already
analyzed epochs from the created epochSequences
---
.../MultiDimensionalRewardUnfolding.cpp | 37 +++---------
.../MultiDimensionalRewardUnfolding.h | 11 ++--
.../helper/rewardbounded/QuantileHelper.cpp | 60 +++++++++----------
.../helper/rewardbounded/QuantileHelper.h | 1 +
4 files changed, 44 insertions(+), 65 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
index fad4feda2..5a1001b3d 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
@@ -64,8 +64,6 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initialize(std::set<storm::expressions::Variable> const& infinityBoundVariables) {
- maxSolutionsStored = 0;
-
STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives.");
std::vector<Epoch> epochSteps;
initializeObjectives(epochSteps, infinityBoundVariables);
@@ -333,41 +331,27 @@ namespace storm {
}
template<typename ValueType, bool SingleObjectiveMode>
- std::vector<typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochComputationOrder(Epoch const& startEpoch) {
+ std::vector<typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochComputationOrder(Epoch const& startEpoch, bool stopAtComputedEpochs) {
// Perform a DFS to find all the reachable epochs
std::vector<Epoch> dfsStack;
std::set<Epoch, std::function<bool(Epoch const&, Epoch const&)>> collectedEpochs(std::bind(&EpochManager::epochClassZigZagOrder, &epochManager, std::placeholders::_1, std::placeholders::_2));
- collectedEpochs.insert(startEpoch);
- dfsStack.push_back(startEpoch);
+ if (!stopAtComputedEpochs || epochSolutions.count(startEpoch) == 0) {
+ collectedEpochs.insert(startEpoch);
+ dfsStack.push_back(startEpoch);
+ }
while (!dfsStack.empty()) {
Epoch currentEpoch = dfsStack.back();
dfsStack.pop_back();
for (auto const& step : possibleEpochSteps) {
Epoch successorEpoch = epochManager.getSuccessorEpoch(currentEpoch, step);
- /*
- for (auto const& e : collectedEpochs) {
- std::cout << "Comparing " << epochManager.toString(e) << " and " << epochManager.toString(successorEpoch) << std::endl;
- if (epochManager.epochClassZigZagOrder(e, successorEpoch)) {
- std::cout << " " << epochManager.toString(e) << " < " << epochManager.toString(successorEpoch) << std::endl;
- }
- if (epochManager.epochClassZigZagOrder(successorEpoch, e)) {
- std::cout << " " << epochManager.toString(e) << " > " << epochManager.toString(successorEpoch) << std::endl;
+ if (!stopAtComputedEpochs || epochSolutions.count(successorEpoch) == 0) {
+ if (collectedEpochs.insert(successorEpoch).second) {
+ dfsStack.push_back(std::move(successorEpoch));
}
}
- */
- if (collectedEpochs.insert(successorEpoch).second) {
- dfsStack.push_back(std::move(successorEpoch));
- }
}
}
- /*
- std::cout << "Resulting order: ";
- for (auto const& e : collectedEpochs) {
- std::cout << epochManager.toString(e) << ", ";
- }
- std::cout << std::endl;
- */
return std::vector<Epoch>(collectedEpochs.begin(), collectedEpochs.end());
}
@@ -629,8 +613,6 @@ namespace storm {
epochModel.objectiveRewardFilter.push_back(storm::utility::vector::filterZero(objRewards));
epochModel.objectiveRewardFilter.back().complement();
}
-
- epochModelSizes.push_back(epochModel.epochMatrix.getRowGroupCount());
}
@@ -846,9 +828,6 @@ namespace storm {
solution.productStateToSolutionVectorMap = productStateToEpochModelInStateMap;
solution.solutions = std::move(inStateSolutions);
epochSolutions[currentEpoch.get()] = std::move(solution);
-
- maxSolutionsStored = std::max((uint64_t) epochSolutions.size(), maxSolutionsStored);
-
}
template<typename ValueType, bool SingleObjectiveMode>
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
index 3844c1af7..db4eb2e6d 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
@@ -54,7 +54,11 @@ namespace storm {
*/
Epoch getStartEpoch(bool setUnknownDimsToBottom = false);
- std::vector<Epoch> getEpochComputationOrder(Epoch const& startEpoch);
+ /*!
+ * Computes a sequence of epochs that need to be analyzed to get a result at the start epoch.
+ * @param stopAtComputedEpochs if set, the search for epochs that need to be computed is stopped at epochs that already have been computed earlier.
+ */
+ std::vector<Epoch> getEpochComputationOrder(Epoch const& startEpoch, bool stopAtComputedEpochs = false);
EpochModel<ValueType, SingleObjectiveMode>& setCurrentEpoch(Epoch const& epoch);
@@ -138,11 +142,6 @@ namespace storm {
std::vector<Dimension<ValueType>> dimensions;
std::vector<storm::storage::BitVector> objectiveDimensions;
-
-
- storm::utility::Stopwatch swInit, swFindSol, swInsertSol, swSetEpoch, swSetEpochClass, swAux1, swAux2, swAux3, swAux4;
- std::vector<uint64_t> epochModelSizes;
- uint64_t maxSolutionsStored;
};
}
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index c6bdda24e..ea7d5ad06 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -186,6 +186,7 @@ namespace storm {
numCheckedEpochs = 0;
numPrecisionRefinements = 0;
swEpochAnalysis.reset();
+ swExploration.reset();
cachedSubQueryResults.clear();
std::vector<std::vector<ValueType>> result;
@@ -217,7 +218,8 @@ namespace storm {
if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl;
std::cout << "Number of required precision refinements: " << numPrecisionRefinements << std::endl;
- std::cout << "Time for epoch model analysis: " << swEpochAnalysis << " seconds." << std::endl;
+ std::cout << "Time for epoch exploration: " << swExploration << " seconds." << std::endl;
+ std::cout << "\tTime for epoch model analysis: " << swEpochAnalysis << " seconds." << std::endl;
}
return result;
}
@@ -359,8 +361,7 @@ namespace storm {
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
std::vector<ValueType> x, b;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
- std::set<EpochManager::Epoch> checkedEpochs;
-
+ swExploration.start();
bool progress = true;
for (CostLimit candidateCostLimitSum(0); progress; ++candidateCostLimitSum.get()) {
CostLimits currentCandidate(satCostLimits.dimension(), CostLimit(0));
@@ -389,42 +390,41 @@ namespace storm {
++costLimitIt;
}
STORM_LOG_DEBUG("Checking start epoch " << rewardUnfolding.getEpochManager().toString(startEpoch) << ".");
- auto epochSequence = rewardUnfolding.getEpochComputationOrder(startEpoch);
+ auto epochSequence = rewardUnfolding.getEpochComputationOrder(startEpoch, true);
for (auto const& epoch : epochSequence) {
- if (checkedEpochs.count(epoch) == 0) {
- ++numCheckedEpochs;
- swEpochAnalysis.start();
- checkedEpochs.insert(epoch);
- auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env,boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
- swEpochAnalysis.stop();
+ ++numCheckedEpochs;
+ swEpochAnalysis.start();
+ auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
+ rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env,boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+ swEpochAnalysis.stop();
- CostLimits epochAsCostLimits;
- if (translateEpochToCostLimits(epoch, startEpoch, consideredDimensions, lowerBoundedDimensions, rewardUnfolding.getEpochManager(), epochAsCostLimits)) {
- ValueType currValue = rewardUnfolding.getInitialStateResult(epoch);
- bool propertySatisfied;
- if (env.solver().isForceSoundness()) {
- ValueType sumOfEpochDimensions = storm::utility::convertNumber<ValueType>(rewardUnfolding.getEpochManager().getSumOfDimensions(epoch) + 1);
- auto lowerUpperValue = getLowerUpperBound(env, sumOfEpochDimensions, currValue);
- propertySatisfied = boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.first);
- if (propertySatisfied != boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.second)) {
- // unclear result due to insufficient precision.
- return false;
- }
- } else {
- propertySatisfied = boundedUntilOperator.getBound().isSatisfied(currValue);
- }
- if (propertySatisfied) {
- satCostLimits.insert(epochAsCostLimits);
- } else {
- unsatCostLimits.insert(epochAsCostLimits);
+ CostLimits epochAsCostLimits;
+ if (translateEpochToCostLimits(epoch, startEpoch, consideredDimensions, lowerBoundedDimensions, rewardUnfolding.getEpochManager(), epochAsCostLimits)) {
+ ValueType currValue = rewardUnfolding.getInitialStateResult(epoch);
+ bool propertySatisfied;
+ if (env.solver().isForceSoundness()) {
+ ValueType sumOfEpochDimensions = storm::utility::convertNumber<ValueType>(rewardUnfolding.getEpochManager().getSumOfDimensions(epoch) + 1);
+ auto lowerUpperValue = getLowerUpperBound(env, sumOfEpochDimensions, currValue);
+ propertySatisfied = boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.first);
+ if (propertySatisfied != boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.second)) {
+ // unclear result due to insufficient precision.
+ swExploration.stop();
+ return false;
}
+ } else {
+ propertySatisfied = boundedUntilOperator.getBound().isSatisfied(currValue);
+ }
+ if (propertySatisfied) {
+ satCostLimits.insert(epochAsCostLimits);
+ } else {
+ unsatCostLimits.insert(epochAsCostLimits);
}
}
}
}
} while (getNextCandidateCostLimit(candidateCostLimitSum, currentCandidate));
}
+ swExploration.stop();
return true;
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
index 63a06e34e..e4377616d 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
@@ -53,6 +53,7 @@ namespace storm {
mutable uint64_t numCheckedEpochs;
mutable uint64_t numPrecisionRefinements;
mutable storm::utility::Stopwatch swEpochAnalysis;
+ mutable storm::utility::Stopwatch swExploration;
};
}
}
From c40ecae2e61d13e169d3f9456b18bc1811a5d902 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Fri, 15 Mar 2019 09:59:29 +0100
Subject: [PATCH 09/11] Implemented quantiles for DTMCs.
---
.../prctl/SparseDtmcPrctlModelChecker.cpp | 32 ++++++++++++++++++-
.../prctl/SparseDtmcPrctlModelChecker.h | 1 +
.../prctl/SparseMdpPrctlModelChecker.cpp | 2 +-
.../helper/rewardbounded/QuantileHelper.cpp | 16 ++++++++--
4 files changed, 47 insertions(+), 4 deletions(-)
diff --git a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index f561d3e98..7c11c31b8 100644
--- a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -7,9 +7,11 @@
#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h"
#include "storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h"
#include "storm/modelchecker/csl/helper/SparseCtmcCslHelper.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h"
#include "storm/logic/FragmentSpecification.h"
@@ -31,7 +33,13 @@ namespace storm {
template<typename SparseDtmcModelType>
bool SparseDtmcPrctlModelChecker<SparseDtmcModelType>::canHandle(CheckTask<storm::logic::Formula, ValueType> const& checkTask) const {
storm::logic::Formula const& formula = checkTask.getFormula();
- return formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setTotalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true).setTimeOperatorsAllowed(true).setReachbilityTimeFormulasAllowed(true));
+ if (formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setTotalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true).setTimeOperatorsAllowed(true).setReachbilityTimeFormulasAllowed(true))) {
+ return true;
+ } else if (formula.isInFragment(storm::logic::quantiles())) {
+ if (this->getModel().getInitialStates().getNumberOfSetBits() > 1) return false;
+ return true;
+ }
+ return false;
}
template<typename SparseDtmcModelType>
@@ -184,6 +192,28 @@ namespace storm {
return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
}
+
+ template<>
+ std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
+ STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Quantiles for parametric models are not implemented.");
+ }
+
+ template<typename SparseDtmcModelType>
+ std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
+ STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for the initial states of a model.");
+ STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, "Quantiles not supported on models with multiple initial states.");
+ uint64_t initialState = *this->getModel().getInitialStates().begin();
+
+ helper::rewardbounded::QuantileHelper<SparseDtmcModelType> qHelper(this->getModel(), checkTask.getFormula());
+ auto res = qHelper.computeQuantile(env);
+
+ if (res.size() == 1 && res.front().size() == 1) {
+ return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(initialState, std::move(res.front().front())));
+ } else {
+ return std::unique_ptr<CheckResult>(new ExplicitParetoCurveCheckResult<ValueType>(initialState, std::move(res)));
+ }
+ }
+
template class SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>>;
#ifdef STORM_HAVE_CARL
diff --git a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index 1aafc0efb..1001b5572 100644
--- a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -34,6 +34,7 @@ namespace storm {
virtual std::unique_ptr<CheckResult> computeConditionalRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::ConditionalFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeReachabilityTimes(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::EventuallyFormula, ValueType> const& checkTask) override;
+ virtual std::unique_ptr<CheckResult> checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) override;
};
} // namespace modelchecker
diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 759d7a976..bd2dfd326 100644
--- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -228,7 +228,7 @@ namespace storm {
template<typename SparseMdpModelType>
std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
- STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for models with a single initial states.");
+ STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for the initial states of a model.");
STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, "Quantiles not supported on models with multiple initial states.");
uint64_t initialState = *this->getModel().getInitialStates().begin();
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index ea7d5ad06..45512006c 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -8,6 +8,7 @@
#include "storm/environment/solver/MinMaxSolverEnvironment.h"
#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/Dtmc.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
#include "storm/storage/expressions/ExpressionManager.h"
#include "storm/storage/expressions/Expressions.h"
@@ -360,7 +361,12 @@ namespace storm {
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
std::vector<ValueType> x, b;
- std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+ std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver; // Needed for MDP
+ std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSolver; // Needed for DTMC
+ if (!model.isNondeterministicModel()) {
+ rewardUnfolding.setEquationSystemFormatForEpochModel(storm::solver::GeneralLinearEquationSolverFactory<ValueType>().getEquationProblemFormat(env));
+ }
+
swExploration.start();
bool progress = true;
for (CostLimit candidateCostLimitSum(0); progress; ++candidateCostLimitSum.get()) {
@@ -395,7 +401,11 @@ namespace storm {
++numCheckedEpochs;
swEpochAnalysis.start();
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
- rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env,boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+ if (model.isNondeterministicModel()) {
+ rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+ } else {
+ rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, x, b, linEqSolver, lowerBound, upperBound));
+ }
swEpochAnalysis.stop();
CostLimits epochAsCostLimits;
@@ -431,6 +441,8 @@ namespace storm {
template class QuantileHelper<storm::models::sparse::Mdp<double>>;
template class QuantileHelper<storm::models::sparse::Mdp<storm::RationalNumber>>;
+ template class QuantileHelper<storm::models::sparse::Dtmc<double>>;
+ template class QuantileHelper<storm::models::sparse::Dtmc<storm::RationalNumber>>;
}
}
From 1ae0200b51432c3ae9598d41579d666d05752fa9 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Fri, 15 Mar 2019 10:00:54 +0100
Subject: [PATCH 10/11] Quantiles: fixed some bugs related to one or three
dimensional quantile queries.
---
.../helper/rewardbounded/CostLimitClosure.cpp | 40 +++++++++++++++++++
.../helper/rewardbounded/CostLimitClosure.h | 6 ++-
.../helper/rewardbounded/QuantileHelper.cpp | 11 +++--
3 files changed, 52 insertions(+), 5 deletions(-)
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
index fefc20299..657703b93 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
@@ -1,6 +1,9 @@
+#include <storm/utility/ExpressionHelper.h>
#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/IllegalArgumentException.h"
+#include "storm/utility/solver.h"
+#include "storm/solver/SmtSolver.h"
namespace storm {
namespace modelchecker {
@@ -140,6 +143,43 @@ namespace storm {
uint64_t CostLimitClosure::dimension() const {
return downwardDimensions.size();
}
+
+ bool CostLimitClosure::unionFull(CostLimitClosure const& first, CostLimitClosure const& second) {
+ assert(first.dimension() == second.dimension());
+ uint64_t dimension = first.dimension();
+ auto manager = std::make_shared<storm::expressions::ExpressionManager>();
+ auto solver = storm::utility::solver::getSmtSolver(*manager);
+
+ std::vector<storm::expressions::Expression> point;
+ storm::expressions::Expression zero = manager->integer(0);
+ for (uint64_t i = 0; i < dimension; ++i) {
+ point.push_back(manager->declareIntegerVariable("x" + std::to_string(i)).getExpression());
+ solver->add(point.back() >= zero);
+ }
+ for (auto const& cl : {first, second}) {
+ for (auto const& q : cl.getGenerator()) {
+ storm::expressions::Expression pointNotDominated;
+ for (uint64_t i = 0; i < point.size(); ++i) {
+ if (!cl.downwardDimensions.get(i) || !q[i].isInfinity()) {
+ assert(!q[i].isInfinity());
+ storm::expressions::Expression qi = manager->integer(q[i].get());
+ storm::expressions::Expression piNotDominated = cl.downwardDimensions.get(i) ? point[i] > qi : point[i] < qi;
+ if (piNotDominated.isInitialized()) {
+ pointNotDominated = pointNotDominated || piNotDominated;
+ } else {
+ pointNotDominated = piNotDominated;
+ }
+ }
+ }
+ if (pointNotDominated.isInitialized()) {
+ solver->add(pointNotDominated);
+ } else {
+ solver->add(manager->boolean(false));
+ }
+ }
+ }
+ return solver->check() == storm::solver::SmtSolver::CheckResult::Unsat;;
+ }
}
}
}
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
index 06abb4722..de6abb4b0 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
@@ -43,7 +43,11 @@ namespace storm {
std::vector<CostLimits> getDominatingCostLimits(CostLimits const& costLimits) const;
GeneratorType const& getGenerator() const;
uint64_t dimension() const;
-
+
+ /*!
+ * Returns true if the union of the two closures is full, i.e., contains every point.
+ */
+ static bool unionFull(CostLimitClosure const& first, CostLimitClosure const& second);
private:
/// The dimensions that are downwards closed, i.e., if x is in the closure, then also all y <= x
diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
index 45512006c..b39eb2b27 100644
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
@@ -273,10 +273,10 @@ namespace storm {
++i;
}
}
- if (subQueryComplement) {
- unsatCostLimits.insert(initPoint);
- } else {
+ if (subQueryComplement == complementaryQuery) {
satCostLimits.insert(initPoint);
+ } else {
+ unsatCostLimits.insert(initPoint);
}
}
}
@@ -298,7 +298,7 @@ namespace storm {
cachedSubQueryResults.emplace(cacheKey, result);
return result;
}
- STORM_LOG_WARN("Restarting quantile computation due to insufficient precision.");
+ STORM_LOG_WARN("Restarting quantile computation after " << swExploration << " seconds due to insufficient precision.");
++numPrecisionRefinements;
increasePrecision(env);
}
@@ -433,6 +433,9 @@ namespace storm {
}
}
} while (getNextCandidateCostLimit(candidateCostLimitSum, currentCandidate));
+ if (!progress) {
+ progress = !CostLimitClosure::unionFull(satCostLimits, unsatCostLimits);
+ }
}
swExploration.stop();
return true;
From d24f61ded6c265ce39e2075cea5d3961a5bca2dc Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Fri, 15 Mar 2019 10:01:29 +0100
Subject: [PATCH 11/11] Added tests for quantiles.
---
.../testfiles/dtmc/quantiles_simple_dtmc.pm | 20 ++
.../testfiles/mdp/quantiles_firewire.nm | 80 +++++
.../testfiles/mdp/quantiles_resources.nm | 72 +++++
.../testfiles/mdp/quantiles_simple_mdp.nm | 31 ++
.../storm/modelchecker/QuantileQueryTest.cpp | 302 ++++++++++++++++++
5 files changed, 505 insertions(+)
create mode 100644 resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm
create mode 100644 resources/examples/testfiles/mdp/quantiles_firewire.nm
create mode 100644 resources/examples/testfiles/mdp/quantiles_resources.nm
create mode 100644 resources/examples/testfiles/mdp/quantiles_simple_mdp.nm
create mode 100644 src/test/storm/modelchecker/QuantileQueryTest.cpp
diff --git a/resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm b/resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm
new file mode 100644
index 000000000..7551aa7a6
--- /dev/null
+++ b/resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm
@@ -0,0 +1,20 @@
+
+dtmc
+
+module main
+ s : [0..3] init 0;
+
+ [] s=0 -> 0.4 : (s'=1) + 0.4 : (s'=2) + 0.2 : (s'=3);
+ [a] s=1 -> 1 : (s'=0);
+ [b] s=2 -> 1 : (s'=0);
+ [] s=3 -> 1 : (s'=3);
+endmodule
+
+rewards "first"
+ [a] true : 1;
+endrewards
+
+rewards "second"
+ [b] true : 2;
+endrewards
+
diff --git a/resources/examples/testfiles/mdp/quantiles_firewire.nm b/resources/examples/testfiles/mdp/quantiles_firewire.nm
new file mode 100644
index 000000000..447f5abf0
--- /dev/null
+++ b/resources/examples/testfiles/mdp/quantiles_firewire.nm
@@ -0,0 +1,80 @@
+// integer semantics version of abstract firewire protocol
+// gxn 23/05/2001
+
+mdp
+
+// wire delay
+const int delay;
+
+// probability of choosing fast and slow
+const double fast = 0.5;
+const double slow = 1-fast;
+
+// largest constant the clock of the system is compared to
+const int kx = 167;
+
+module abstract_firewire
+
+ // clock
+ x : [0..kx+1];
+
+ // local state
+ s : [0..9];
+ // 0 -start_start
+ // 1 -fast_start
+ // 2 -start_fast
+ // 3 -start_slow
+ // 4 -slow_start
+ // 5 -fast_fast
+ // 6 -fast_slow
+ // 7 -slow_fast
+ // 8 -slow_slow
+ // 9 -done
+
+ // initial state
+ [time] s=0 & x<delay -> (x'=min(x+1,kx+1));
+ [round] s=0 -> fast : (s'=1) + slow : (s'=4);
+ [round] s=0 -> fast : (s'=2) + slow : (s'=3);
+ // fast_start
+ [time] s=1 & x<delay -> (x'=min(x+1,kx+1));
+ [] s=1 -> fast : (s'=5) & (x'=0) + slow : (s'=6) & (x'=0);
+ // start_fast
+ [time] s=2 & x<delay -> (x'=min(x+1,kx+1));
+ [] s=2 -> fast : (s'=5) & (x'=0) + slow : (s'=7) & (x'=0);
+ // start_slow
+ [time] s=3 & x<delay -> (x'=min(x+1,kx+1));
+ [] s=3 -> fast : (s'=6) & (x'=0) + slow : (s'=8) & (x'=0);
+ // slow_start
+ [time] s=4 & x<delay -> (x'=min(x+1,kx+1));
+ [] s=4 -> fast : (s'=7) & (x'=0) + slow : (s'=8) & (x'=0);
+ // fast_fast
+ [time] s=5 & (x<85) -> (x'=min(x+1,kx+1));
+ [] s=5 & (x>=76) -> (s'=0) & (x'=0);
+ [] s=5 & (x>=76-delay) -> (s'=9) & (x'=0);
+ // fast_slow
+ [time] s=6 & x<167 -> (x'=min(x+1,kx+1));
+ [] s=6 & x>=159-delay -> (s'=9) & (x'=0);
+ // slow_fast
+ [time] s=7 & x<167 -> (x'=min(x+1,kx+1));
+ [] s=7 & x>=159-delay -> (s'=9) & (x'=0);
+ // slow_slow
+ [time] s=8 & x<167 -> (x'=min(x+1,kx+1));
+ [] s=8 & x>=159 -> (s'=0) & (x'=0);
+ [] s=8 & x>=159-delay -> (s'=9) & (x'=0);
+ // done
+ [] s=9 -> (s'=s);
+
+endmodule
+
+// labels
+label "done" = s=9;
+
+//reward structures
+// time
+rewards "time"
+ [time] true : 1;
+endrewards
+// number of rounds
+rewards "rounds"
+ [round] true : 1;
+endrewards
diff --git a/resources/examples/testfiles/mdp/quantiles_resources.nm b/resources/examples/testfiles/mdp/quantiles_resources.nm
new file mode 100644
index 000000000..6061465e4
--- /dev/null
+++ b/resources/examples/testfiles/mdp/quantiles_resources.nm
@@ -0,0 +1,72 @@
+mdp
+
+const int WIDTH = 5;
+const int HEIGHT = 5;
+const int XINIT = 3;
+const int YINIT = 1;
+
+const double pAttack = 1/10;
+
+formula left_of_gold = x=2 & y=5;
+formula right_of_gold = x=4 & y=5;
+formula below_of_gold = (x=3 & y=4);
+formula above_of_gold = false;
+formula left_of_gem = (x=4 & y=4);
+formula right_of_gem = false;
+formula below_of_gem = (x=5 & y=3);
+formula above_of_gem = (x=5 & y=5);
+formula left_of_home = x=2 & y=1;
+formula right_of_home = x=4 & y=1;
+formula above_of_home = x=3 & y=2;
+formula below_of_home = false;
+formula left_of_enemy = (x=3 & y=5) | (x=2 & y=4);
+formula right_of_enemy = (x=5 & y=5) | (x=4 & y=4);
+formula above_of_enemy = x=3 & y=5;
+formula below_of_enemy = (x=3 & y=3) | (x=4 & y=4);
+
+module robot
+
+ gold : bool init false;
+ gem : bool init false;
+ attacked : bool init false;
+
+ x : [1..WIDTH] init XINIT;
+ y : [1..HEIGHT] init YINIT;
+
+ [right] !left_of_enemy & x<WIDTH -> (attacked'=false) & (x'=x+1) & (gold' = (gold & !left_of_home) | left_of_gold) & (gem' = (gem & !left_of_home) | left_of_gem);
+ [left] !right_of_enemy & x>1 -> (attacked'=false) & (x'=x-1) & (gold' = (gold & !right_of_home) | right_of_gold) & (gem' = (gem & !right_of_home) | right_of_gem);
+ [top] !below_of_enemy & y<HEIGHT -> (attacked'=false) & (y'=y+1) & (gold' = (gold & !below_of_home) | below_of_gold) & (gem' = (gem & !below_of_home) | below_of_gem);
+ [down] !above_of_enemy & y>1 -> (attacked'=false) & (y'=y-1) & (gold' = (gold & !above_of_home) | above_of_gold) & (gem' = (gem & !above_of_home) | above_of_gem);
+
+ [right] left_of_enemy & x<WIDTH -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (x'=x+1) & (gold' = (gold & !left_of_home) | left_of_gold) & (gem' = (gem & !left_of_home) | left_of_gem);
+ [left] right_of_enemy & x>1 -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (x'=x-1) & (gold' = (gold & !right_of_home) | right_of_gold) & (gem' = (gem & !right_of_home) | right_of_gem);
+ [top] below_of_enemy & y<HEIGHT -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (y'=y+1) & (gold' = (gold & !below_of_home) | below_of_gold) & (gem' = (gem & !below_of_home) | below_of_gem);
+ [down] above_of_enemy & y>1 -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (y'=y-1) & (gold' = (gold & !above_of_home) | above_of_gold) & (gem' = (gem & !above_of_home) | above_of_gem);
+endmodule
+
+rewards "attacks"
+ attacked : 1;
+endrewards
+
+rewards "steps"
+ [left] true : 1;
+ [right] true : 1;
+ [top] true : 1;
+ [down] true : 1;
+endrewards
+
+rewards "gold"
+ [right] left_of_home & gold : 1;
+ [left] right_of_home & gold : 1;
+ [top] below_of_home & gold : 1;
+ [down] above_of_home & gold : 1;
+endrewards
+
+rewards "gem"
+ [right] left_of_home & gem : 1;
+ [left] right_of_home & gem : 1;
+ [top] below_of_home & gem : 1;
+ [down] above_of_home & gem: 1;
+endrewards
+
+
\ No newline at end of file
diff --git a/resources/examples/testfiles/mdp/quantiles_simple_mdp.nm b/resources/examples/testfiles/mdp/quantiles_simple_mdp.nm
new file mode 100644
index 000000000..011bc6232
--- /dev/null
+++ b/resources/examples/testfiles/mdp/quantiles_simple_mdp.nm
@@ -0,0 +1,31 @@
+
+mdp
+
+module main
+ s : [0..6] init 6;
+
+ [] s=6 -> 1 : (s'=0);
+ [] s=6 -> 1 : (s'=5);
+ [a] s=5 -> 1 : (s'=5);
+ [b] s=5 -> 1 : (s'=1);
+ [] s=0 -> 1/10: (s'=1) + 9/10: (s'=3);
+ [] s=0 -> 1/10: (s'=1) + 9/10: (s'=4);
+ [a] s=3 -> 1: (s'=0);
+ [b] s=3 -> 1: (s'=0);
+ [] s=1 -> 1: (s'=2);
+ [] s=2 -> 1: (s'=2);
+endmodule
+
+rewards "first"
+ [a] true : 1;
+endrewards
+
+rewards "second"
+ [b] s=3 : 1;
+ [b] s=5 : 2;
+endrewards
+
+rewards "third"
+ true : 0.7;
+endrewards
+
diff --git a/src/test/storm/modelchecker/QuantileQueryTest.cpp b/src/test/storm/modelchecker/QuantileQueryTest.cpp
new file mode 100644
index 000000000..c6c718856
--- /dev/null
+++ b/src/test/storm/modelchecker/QuantileQueryTest.cpp
@@ -0,0 +1,302 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "test/storm_gtest.h"
+
+#include "storm/api/builder.h"
+#include "storm-parsers/api/model_descriptions.h"
+#include "storm/api/properties.h"
+#include "storm-parsers/api/properties.h"
+#include "storm/parser/CSVParser.h"
+
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
+#include "storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+#include "storm/environment/solver/MinMaxSolverEnvironment.h"
+#include "storm/logic/Formulas.h"
+#include "storm/storage/jani/Property.h"
+
+namespace {
+
+ enum class MdpEngine {PrismSparse, JaniSparse, JitSparse, Hybrid, PrismDd, JaniDd};
+
+ class UnsoundEnvironment {
+ public:
+ typedef double ValueType;
+ static storm::Environment createEnvironment() {
+ storm::Environment env;
+ env.solver().minMax().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
+ return env;
+ }
+ };
+
+ class SoundEnvironment {
+ public:
+ typedef double ValueType;
+ static storm::Environment createEnvironment() {
+ storm::Environment env;
+ env.solver().setForceSoundness(true);
+ return env;
+ }
+ };
+
+ class ExactEnvironment {
+ public:
+ typedef storm::RationalNumber ValueType;
+ static storm::Environment createEnvironment() {
+ storm::Environment env;
+ return env;
+ }
+ };
+
+ template<typename TestType>
+ class QuantileQueryTest : public ::testing::Test {
+ public:
+ typedef typename TestType::ValueType ValueType;
+ QuantileQueryTest() : _environment(TestType::createEnvironment()) {}
+ storm::Environment const& env() const { return _environment; }
+ ValueType parseNumber(std::string const& input) const {
+ if (input.find("inf") != std::string::npos) {
+ return storm::utility::infinity<ValueType>();
+ }
+ return storm::utility::convertNumber<ValueType>(input);
+ }
+
+ template <typename MT>
+ std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> buildModelFormulas(std::string const& pathToPrismFile, std::string const& formulasAsString, std::string const& constantDefinitionString = "") const {
+ std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> result;
+ storm::prism::Program program = storm::api::parseProgram(pathToPrismFile);
+ program = storm::utility::prism::preprocess(program, constantDefinitionString);
+ result.second = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program));
+ result.first = storm::api::buildSparseModel<ValueType>(program, result.second)->template as<MT>();
+ return result;
+ }
+
+ std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> getTasks(std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas) const {
+ std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> result;
+ for (auto const& f : formulas) {
+ result.emplace_back(*f, true); // only initial states are relevant
+ }
+ return result;
+ }
+
+ template <typename MT>
+ typename std::enable_if<std::is_same<MT, storm::models::sparse::Mdp<ValueType>>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
+ createModelChecker(std::shared_ptr<MT> const& model) const {
+ return std::make_shared<storm::modelchecker::SparseMdpPrctlModelChecker<MT>>(*model);
+ }
+ template <typename MT>
+ typename std::enable_if<std::is_same<MT, storm::models::sparse::Dtmc<ValueType>>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
+ createModelChecker(std::shared_ptr<MT> const& model) const {
+ return std::make_shared<storm::modelchecker::SparseDtmcPrctlModelChecker<MT>>(*model);
+ }
+
+ std::pair<bool, std::string> compareResult(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, std::unique_ptr<storm::modelchecker::CheckResult>& result, std::vector<std::string> const& expected) {
+ bool equal = true;
+ std::string errorMessage = "";
+ auto filter = getInitialStateFilter(model);
+ result->filter(*filter);
+ std::vector<std::vector<ValueType>> resultPoints;
+ if (result->isExplicitParetoCurveCheckResult()) {
+ resultPoints = result->asExplicitParetoCurveCheckResult<ValueType>().getPoints();
+ } else {
+ if (!result->isExplicitQuantitativeCheckResult()) {
+ return { false, "The CheckResult has unexpected type."};
+ }
+ resultPoints = {{ result->asExplicitQuantitativeCheckResult<ValueType>().getMax() }};
+ }
+ std::vector<std::vector<ValueType>> expectedPoints;
+ for (auto const& pointAsString : expected) {
+ std::vector<ValueType> point;
+ for (auto const& entry : storm::parser::parseCommaSeperatedValues(pointAsString)) {
+ point.push_back(parseNumber(entry));
+ }
+ expectedPoints.push_back(std::move(point));
+ }
+ for (auto const& resPoint : resultPoints) {
+ bool contained = false;
+ for (auto const& expPoint : expectedPoints) {
+ if (resPoint == expPoint) {
+ contained = true;
+ break;
+ }
+ }
+ if (!contained) {
+ equal = false;
+ errorMessage += "The result " + storm::utility::vector::toString(resPoint) + " is not contained in the expected solution. ";
+ }
+ }
+ for (auto const& expPoint : expectedPoints) {
+ bool contained = false;
+ for (auto const& resPoint : resultPoints) {
+ if (resPoint == expPoint) {
+ contained = true;
+ break;
+ }
+ }
+ if (!contained) {
+ equal = false;
+ errorMessage += "The expected " + storm::utility::vector::toString(expPoint) + " is not contained in the obtained solution. ";
+ }
+ }
+ return {equal, errorMessage};
+ }
+
+ private:
+ storm::Environment _environment;
+
+ std::unique_ptr<storm::modelchecker::QualitativeCheckResult> getInitialStateFilter(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model) const {
+ return std::make_unique<storm::modelchecker::ExplicitQualitativeCheckResult>(model->getInitialStates());
+ }
+ };
+
+ typedef ::testing::Types<
+ UnsoundEnvironment,
+ SoundEnvironment,
+ ExactEnvironment
+ > TestingTypes;
+
+ TYPED_TEST_CASE(QuantileQueryTest, TestingTypes);
+
+
+ TYPED_TEST(QuantileQueryTest, simple_Dtmc) {
+ typedef storm::models::sparse::Dtmc<typename TestFixture::ValueType> ModelType;
+
+ std::string formulasString = "quantile(max A, max B, P>0.95 [F{\"first\"}<=A,{\"second\"}<=B s=3]);\n";
+
+ auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/dtmc/quantiles_simple_dtmc.pm", formulasString);
+ auto model = std::move(modelFormulas.first);
+ auto tasks = this->getTasks(modelFormulas.second);
+ EXPECT_EQ(4ul, model->getNumberOfStates());
+ EXPECT_EQ(6ul, model->getNumberOfTransitions());
+ auto checker = this->template createModelChecker<ModelType>(model);
+ std::unique_ptr<storm::modelchecker::CheckResult> result;
+ std::vector<std::string> expectedResult;
+ std::pair<bool, std::string> compare;
+ uint64_t taskId = 0;
+
+ expectedResult.clear();
+ expectedResult.push_back("7, 18");
+ expectedResult.push_back("8, 16");
+ expectedResult.push_back("9, 14");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+ }
+
+ TYPED_TEST(QuantileQueryTest, simple_Mdp) {
+ typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
+
+ std::string formulasString = "quantile(B1, B2, Pmax>0.5 [F{\"first\"}>=B1,{\"second\"}>=B2 s=1]);\n";
+ formulasString += "quantile(B1, B2, Pmax>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2 s=1]);\n";
+ formulasString += "quantile(B1, B2, Pmin>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2 s=1]);\n";
+ formulasString += "quantile(B1, Pmax>0.5 [F{\"second\"}>=B1 s=1]);\n";
+ formulasString += "quantile(B1, B2, B3, Pmax>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2,{\"third\"}<=B3 s=1]);\n";
+
+ auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_simple_mdp.nm", formulasString);
+ auto model = std::move(modelFormulas.first);
+ auto tasks = this->getTasks(modelFormulas.second);
+ EXPECT_EQ(7ul, model->getNumberOfStates());
+ EXPECT_EQ(13ul, model->getNumberOfTransitions());
+ auto checker = this->template createModelChecker<ModelType>(model);
+ std::unique_ptr<storm::modelchecker::CheckResult> result;
+ std::vector<std::string> expectedResult;
+ std::pair<bool, std::string> compare;
+ uint64_t taskId = 0;
+
+ expectedResult.clear();
+ expectedResult.push_back(" 0, 6");
+ expectedResult.push_back(" 1, 5");
+ expectedResult.push_back(" 2, 4");
+ expectedResult.push_back(" 3, 3");
+ expectedResult.push_back("inf, 2");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ expectedResult.clear();
+ expectedResult.push_back(" 0, 2");
+ expectedResult.push_back(" 5, 1");
+ expectedResult.push_back(" 6, 0");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ expectedResult.clear();
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ expectedResult.clear();
+ expectedResult.push_back(" 6");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ expectedResult.clear();
+ expectedResult.push_back(" 0, 2, 1.4");
+ expectedResult.push_back(" 5, 1, 9.8");
+ expectedResult.push_back(" 6, 0, 9.8");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+ }
+
+ TYPED_TEST(QuantileQueryTest, firewire) {
+ typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
+
+ std::string formulasString = "quantile(min TIME, min ROUNDS, Pmax>0.95 [F{\"time\"}<=TIME,{\"rounds\"}<=ROUNDS \"done\"]);\n";
+
+ auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_firewire.nm", formulasString, "delay=36");
+ auto model = std::move(modelFormulas.first);
+ auto tasks = this->getTasks(modelFormulas.second);
+ EXPECT_EQ(776ul, model->getNumberOfStates());
+ EXPECT_EQ(1411ul, model->getNumberOfTransitions());
+ auto checker = this->template createModelChecker<ModelType>(model);
+ std::unique_ptr<storm::modelchecker::CheckResult> result;
+ std::vector<std::string> expectedResult;
+ std::pair<bool, std::string> compare;
+ uint64_t taskId = 0;
+
+ expectedResult.clear();
+ expectedResult.push_back("123, 1");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ }
+
+ TYPED_TEST(QuantileQueryTest, resources) {
+ typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
+
+ std::string formulasString = "quantile(max GOLD, max GEM, Pmax>0.95 [F{\"gold\"}>=GOLD,{\"gem\"}>=GEM,{\"steps\"}<=100 true]);\n";
+
+ auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_resources.nm", formulasString);
+ auto model = std::move(modelFormulas.first);
+ auto tasks = this->getTasks(modelFormulas.second);
+ EXPECT_EQ(94ul, model->getNumberOfStates());
+ EXPECT_EQ(326ul, model->getNumberOfTransitions());
+ auto checker = this->template createModelChecker<ModelType>(model);
+ std::unique_ptr<storm::modelchecker::CheckResult> result;
+ std::vector<std::string> expectedResult;
+ std::pair<bool, std::string> compare;
+ uint64_t taskId = 0;
+
+ expectedResult.clear();
+ expectedResult.push_back("0, 10");
+ expectedResult.push_back("1, 9");
+ expectedResult.push_back("4, 8");
+ expectedResult.push_back("7, 7");
+ expectedResult.push_back("8, 4");
+ expectedResult.push_back("9, 2");
+ expectedResult.push_back("10, 0");
+ result = checker->check(this->env(), tasks[taskId++]);
+ compare = this->compareResult(model, result, expectedResult);
+ EXPECT_TRUE(compare.first) << compare.second;
+
+ }
+}