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 CostLimitClosure::getDominatingCostLimits(CostLimits const& costLimits) const { std::vector 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 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 EpochModel& MultiDimensionalRewardUnfolding::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 + template::type> + void MultiDimensionalRewardUnfolding::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 + template::type> + void MultiDimensionalRewardUnfolding::setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const { + STORM_LOG_ASSERT(objIndex < solution.size(), "Invalid objective index " << objIndex << "."); + solution[objIndex] = value; + } + template template::type> std::string MultiDimensionalRewardUnfolding::solutionToString(SolutionType const& solution) const { @@ -858,21 +872,39 @@ namespace storm { template typename MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::getInitialStateResult(Epoch const& epoch) { + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) { + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::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()); + } + } + } + return result; } template @@ -885,11 +917,6 @@ namespace storm { return dimensions.at(dim); } - template - storm::storage::BitVector const& MultiDimensionalRewardUnfolding::getProb1Objectives() const { - return productModel->getProb1Objectives(); - } - template class MultiDimensionalRewardUnfolding; template class MultiDimensionalRewardUnfolding; template class MultiDimensionalRewardUnfolding; 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 getLowerObjectiveBound(uint64_t objectiveIndex = 0); void setSolutionForCurrentEpoch(std::vector&& 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 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::type = 0> void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const; + + template::type = 0> + void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const; + template::type = 0> + void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const; + + + template::type = 0> std::string solutionToString(SolutionType const& solution) const; template::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 - ProductModel::ProductModel(storm::models::sparse::Model const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()), prob1Objectives(objectives.size(), false) { + ProductModel::ProductModel(storm::models::sparse::Model const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector 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 - uint64_t ProductModel::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const { + uint64_t ProductModel::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 @@ -649,8 +648,8 @@ namespace storm { } template - storm::storage::BitVector const& ProductModel::getProb1Objectives() { - return prob1Objectives; + boost::optional const& ProductModel::getProb1InitialStates(uint64_t objectiveIndex) const { + return prob1InitialStates[objectiveIndex]; } template class ProductModel; 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 const& getProb1InitialStates(uint64_t objectiveIndex) const; private: @@ -77,7 +77,7 @@ namespace storm { std::vector productToModelStateMap; std::vector productToMemoryStateMap; std::vector choiceToStateMap; - storm::storage::BitVector prob1Objectives; /// Objectives that are already satisfied in the initial state + std::vector> 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 QuantileHelper::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) { - // Intentionally left empty + // Do all kinds of sanity check. + std::set 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 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 transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector const& transformations) { + std::shared_ptr transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector 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> leftSubformulas, rightSubformulas; @@ -107,7 +124,11 @@ namespace storm { } else { newBoundedUntil = std::make_shared(origBoundedUntil.getLeftSubformula().asSharedPointer(), origBoundedUntil.getRightSubformula().asSharedPointer(), lowerBounds, upperBounds, timeBoundReferences); } - return std::make_shared(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(newBoundedUntil, newOpInfo); } /// Increases the precision of solver results @@ -174,55 +195,26 @@ namespace storm { } template - storm::storage::BitVector QuantileHelper::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 - std::vector> QuantileHelper::computeQuantile(Environment const& env) const { + std::vector> QuantileHelper::computeQuantile(Environment const& env) { + numCheckedEpochs = 0; + numPrecisionRefinements = 0; + cachedSubQueryResults.clear(); + std::vector> 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(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() }}; + // 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 resultPoint(costLimits.size()); + storm::utility::vector::applyPointwise(costLimits, internalResult.second, resultPoint, [](CostLimit const& costLimit, ValueType const& factor) -> ValueType { + if (costLimit.isInfinity()) { + return storm::utility::infinity(); } else { - result = {{ storm::utility::infinity()}}; - } - } 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(costLimit.get()) * factor; + }}); + result.push_back(resultPoint); } if (storm::settings::getModule().isShowStatisticsSet()) { std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl; @@ -231,20 +223,26 @@ 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::pair::ValueType>> QuantileHelper::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(getDimension(), BoundTransformation::None), complementaryQuery); std::set 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 rewardUnfolding(model, boundedUntilOp, infinityVariables); - if (computeQuantile(env, consideredDimensions, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) { + if (computeQuantile(env, consideredDimensions, *boundedUntilOp, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) { std::vector scalingFactors; for (auto const& dim : consideredDimensions) { scalingFactors.push_back(rewardUnfolding.getDimension(dim).scalingFactor); } - return {satCostLimits, scalingFactors}; + std::pair> 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 - bool QuantileHelper::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding& rewardUnfolding) const { + bool QuantileHelper::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding& rewardUnfolding) { auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); auto upperBound = rewardUnfolding.getUpperObjectiveBound(); @@ -341,12 +366,14 @@ namespace storm { std::unique_ptr> minMaxSolver; std::set 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 - void filterDominatedPoints(std::vector>& points, std::vector const& dirs) { - std::vector> 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 - std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment& env) const { - std::vector> result; - - auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); - auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula(); - std::vector dimensionsAsBitVector; - std::vector dimensions; - std::vector optimizationDirections; - std::vector 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> 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()); - } 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 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(), storm::utility::zero()}}; - } else { - result = {{ storm::utility::infinity(), storm::utility::infinity()}}; - } - } 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 - bool QuantileHelper::exploreTwoDimensionalQuantile(Environment const& env, std::vector> const& startEpochValues, std::vector& currentEpochValues, std::vector>& resultPoints) const { - // Initialize some data for easy access - auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); - auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula(); - std::vector dimensionsAsBitVector; - std::vector dimensions; - std::vector optimizationDirections; - std::vector 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 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]); - 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 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 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(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 newResultPoint(2); - newResultPoint[i] = storm::utility::convertNumber(startEpochValues[i].first - 1) * startEpochValues[i].second; - newResultPoint[1-i] = storm::utility::infinity(); - resultPoints.push_back(newResultPoint); - } - } - - filterDominatedPoints(resultPoints, optimizationDirections); - return true; - } - template - boost::optional> QuantileHelper::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 infinityVariables; - for (auto const& d : getOpenDimensions()) { - if (d != dimension) { - infinityVariables.insert(getVariableForDimension(d)); - } - } - auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector(getDimension(), BoundTransformation::None)); - MultiDimensionalRewardUnfolding 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 x, b; - std::unique_ptr> minMaxSolver; - - uint64_t epochValue = 0; - std::set 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 - bool QuantileHelper::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::ValueType QuantileHelper::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 bts(getDimension(), BoundTransformation::None); - std::set 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 rewardUnfolding(model, transformedFormula, infinityVariables); - if (!rewardUnfolding.getProb1Objectives().empty()) { - assert(rewardUnfolding.getProb1Objectives().size() == 1); - // The probability is one. - return storm::utility::one(); - } - // 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 x, b; - std::unique_ptr> 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>; template class QuantileHelper>; 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 +#include + #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> computeQuantile(Environment const& env) const; + std::vector> computeQuantile(Environment const& env); 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::pair> 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& rewardUnfolding); std::vector> computeTwoDimensionalQuantile(Environment& env) const; @@ -66,6 +66,7 @@ namespace storm { ModelType const& model; storm::logic::QuantileFormula const& quantileFormula; + std::map>> cachedSubQueryResults; /// Statistics mutable uint64_t numCheckedEpochs;