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::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::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 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 CostLimitClosure::getDominatingCostLimits(CostLimits const& costLimits) const { + std::vector 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 +#include +#include + +#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 CostLimits; + + class CostLimitClosure { + public: + + struct CostLimitsCompare { + bool operator() (CostLimits const& lhs, CostLimits const& rhs) const; + }; + typedef std::set 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 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 + std::pair::ValueType>> QuantileHelper::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions) const { + // Todo: ask the cache first + + auto boundedUntilOp = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector(getDimension(), BoundTransformation::None)); + std::set 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 rewardUnfolding(model, boundedUntilOp, infinityVariables); + if (computeQuantile(env, consideredDimensions, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) { + std::vector 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(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 + bool QuantileHelper::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding& rewardUnfolding) const { + + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + std::vector x, b; + std::unique_ptr> minMaxSolver; + std::set 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 void filterDominatedPoints(std::vector>& points, std::vector const& dirs) { std::vector> result; @@ -351,13 +531,13 @@ namespace storm { MultiDimensionalRewardUnfolding rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector(getDimension(), BoundTransformation::None))); - + // initialize data that will be needed for each epoch auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); auto upperBound = rewardUnfolding.getUpperObjectiveBound(); std::vector x, b; std::unique_ptr> 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> computeQuantile(Environment const& env) const; private: + + std::pair> 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& rewardUnfolding); + + std::vector> computeTwoDimensionalQuantile(Environment& env) const; bool exploreTwoDimensionalQuantile(Environment const& env, std::vector> const& startEpochValues, std::vector& currentEpochValues, std::vector>& resultPoints) const;