Browse Source

rewardbounded: Various fixes.

tempestpy_adaptions
TimQu 6 years ago
parent
commit
fb7078770d
  1. 26
      src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
  2. 12
      src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
  3. 8
      src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
  4. 2
      src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
  5. 317
      src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
  6. 19
      src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h

26
src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp

@ -278,7 +278,7 @@ namespace storm {
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(upperBoundedDimensions.empty() || !probabilityOperatorFormula.getSubformula().asBoundedUntilFormula().isMultiDimensional(), 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_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);
if (!upperBoundedDimensions.empty()) {
// To not invalidate upper-bounded dimensions, one needs to consider MECS where no reward for such a dimension is collected.
@ -371,6 +371,7 @@ 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())) {
@ -631,19 +632,6 @@ namespace storm {
STORM_LOG_ASSERT(model.isOfType(storm::models::ModelType::Dtmc), "Trying to set the equation problem format although the model is not deterministic.");
epochModel.equationSolverProblemFormat = eqSysFormat;
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<SO, int>::type>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getOneSolution() const {
return storm::utility::one<ValueType>;
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<!SO, int>::type>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getOneSolution() const {
return SolutionType(objectives.size(), storm::utility::one<ValueType>());
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<SO, int>::type>
@ -865,10 +853,7 @@ namespace storm {
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.");
if (!productModel->getProb1Objectives().empty()) {
STORM_LOG_THROW(SingleObjectiveMode, storm::exceptions::NotSupportedException, "One of the objectives is already satisfied in the initial state. This is not implemented in multi-objective mode.");
return getOneSolution();
}
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'
return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]];
}
@ -894,6 +879,11 @@ 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>;

12
src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h

@ -77,6 +77,12 @@ namespace storm {
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:
@ -89,12 +95,6 @@ namespace storm {
void initializeMemoryProduct(std::vector<Epoch> const& epochSteps);
// Returns a solution only consisting of one
template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
SolutionType getOneSolution() const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
SolutionType getOneSolution() const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>

8
src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp

@ -92,9 +92,8 @@ namespace storm {
computeReachableStatesInEpochClasses();
}
template<typename ValueType>
storm::storage::MemoryStructure ProductModel<ValueType>::computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const {
storm::storage::MemoryStructure ProductModel<ValueType>::computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) {
storm::modelchecker::SparsePropositionalModelChecker<storm::models::sparse::Model<ValueType>> mc(model);
@ -649,6 +648,11 @@ namespace storm {
return getProductState(getModelState(productState), transformMemoryState(getMemoryState(productState), epochClass, predecessorMemoryState));
}
template<typename ValueType>
storm::storage::BitVector const& ProductModel<ValueType>::getProb1Objectives() {
return prob1Objectives;
}
template class ProductModel<double>;
template class ProductModel<storm::RationalNumber>;
}

2
src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h

@ -52,7 +52,7 @@ namespace storm {
private:
storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const;
storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
std::vector<MemoryState> computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const;

317
src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp

@ -5,6 +5,8 @@
#include <memory>
#include <boost/optional.hpp>
#include "storm/environment/solver/MinMaxSolverEnvironment.h"
#include "storm/models/sparse/Mdp.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
#include "storm/storage/expressions/ExpressionManager.h"
@ -104,6 +106,33 @@ namespace storm {
return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, boundedUntilOperator.getOperatorInformation());
}
/// Increases the precision of solver results
void increasePrecision(storm::Environment& env) {
STORM_LOG_DEBUG("Increasing precision of underlying solver.");
auto factor = storm::utility::convertNumber<storm::RationalNumber, std::string>("0.1");
env.solver().setLinearEquationSolverPrecision(static_cast<storm::RationalNumber>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get() * factor));
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
template<typename ValueType>
std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& value, bool minMax = true) {
ValueType prec;
bool relative;
if (minMax) {
prec = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision());
relative = env.solver().minMax().getRelativeTerminationCriterion();
} else {
prec = storm::utility::convertNumber<ValueType>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get());
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)));
} else {
return std::make_pair<ValueType, ValueType>(value - prec, value + prec);
}
}
template<typename ModelType>
uint64_t QuantileHelper<ModelType>::getDimension() const {
return quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().getDimension();
@ -158,15 +187,21 @@ namespace storm {
template<typename ModelType>
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) {
std::vector<std::vector<ValueType>> result;
Environment envCpy = env; // It might be necessary to increase the precision during the computation
if (getOpenDimensions().getNumberOfSetBits() == 1) {
uint64_t dimension = *getOpenDimensions().begin();
auto const& boundedUntilOperator = quantileFormula.getSubformula().asProbabilityOperatorFormula();
bool zeroSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeLimitValue(env, storm::storage::BitVector(getDimension(), false)));
bool infSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeLimitValue(env, getOpenDimensions()));
bool zeroSatisfiesFormula = checkLimitValue(envCpy, storm::storage::BitVector(getDimension(), false));
bool infSatisfiesFormula = checkLimitValue(envCpy, getOpenDimensions());
if (zeroSatisfiesFormula != infSatisfiesFormula) {
auto quantileRes = computeQuantileForDimension(env, dimension);
result = {{storm::utility::convertNumber<ValueType>(quantileRes.first) * quantileRes.second}};
while (true) {
auto quantileRes = computeQuantileForDimension(envCpy, dimension);
if (quantileRes) {
result = {{storm::utility::convertNumber<ValueType>(quantileRes->first) * quantileRes->second}};
break;
}
increasePrecision(envCpy);
}
} else if (zeroSatisfiesFormula) {
// thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula
if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) {
@ -179,7 +214,7 @@ namespace storm {
result = {{}};
}
} else if (getOpenDimensions().getNumberOfSetBits() == 2) {
result = computeTwoDimensionalQuantile(env);
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.");
}
@ -218,7 +253,7 @@ namespace storm {
}
template<typename ModelType>
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment const& env) {
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment& env) {
std::vector<std::vector<ValueType>> result;
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
@ -240,15 +275,13 @@ namespace storm {
if (lowerCostBounds[0] == lowerCostBounds[1]) {
// TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0])
bool infInfProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeLimitValue(env, storm::storage::BitVector(getDimension(), false)));
bool zeroZeroProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeLimitValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1]));
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;
std::vector<std::vector<int64_t>> resultPoints;
const int64_t infinity = std::numeric_limits<int64_t>().max(); // use this value to represent infinity in a result point
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 = probOpFormula.getBound().isSatisfied(computeLimitValue(env, dimensionsAsBitVector[1-i]));
bool zeroInfSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[1-i]);
std::cout << "Formula sat is " << zeroInfSatisfiesFormula << " and lower bound is " << storm::logic::isLowerBound(probabilityBound) << std::endl;
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
@ -256,7 +289,7 @@ namespace storm {
} else {
// Compute quantile where 1-i is set to infinity
std::cout << "Computing quantile for single dimension " << dimensions[i] << std::endl;
singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]));
singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]).get());
std::cout << ".. Result is " << singleQuantileValues.back().first << std::endl;
if (!storm::solver::minimize(optimizationDirections[i])) {
// When maximizing bounds, the computed single dimensional quantile value is sat for all values of the other bound.
@ -270,93 +303,10 @@ namespace storm {
}
}
std::cout << "Found starting point. Beginning 2D exploration" << std::endl;
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;
std::vector<int64_t> epochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0).
std::cout << "Starting quantile exploration with: (" << epochValues[0] << ", " << epochValues[1] << ")." << std::endl;
std::set<typename EpochManager::Epoch> exploredEpochs;
while (true) {
auto currEpoch = rewardUnfolding.getStartEpoch(true);
for (uint64_t i = 0; i < 2; ++i) {
if (epochValues[i] >= 0) {
rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) epochValues[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));
exploredEpochs.insert(epoch);
}
}
ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
std::cout << "Numeric result for epoch " << rewardUnfolding.getEpochManager().toString(currEpoch) << " is " << currValue << std::endl;
bool 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!
resultPoints.push_back(epochValues);
std::cout << "Found another result point: (" << epochValues[0] << ", " << epochValues[1] << ")." << std::endl;
if (epochValues[1] == singleQuantileValues[1].first) {
break;
} else {
++epochValues[0];
epochValues[1] = singleQuantileValues[1].first;
}
} else {
++epochValues[1];
}
}
// Translate the result points to the 'original' domain
for (auto& p : resultPoints) {
std::vector<ValueType> convertedP;
for (uint64_t i = 0; i < 2; ++i) {
if (p[i] == infinity) {
convertedP.push_back(storm::utility::infinity<ValueType>());
} else {
if (lowerCostBounds[i]) {
// Translate > to >=
++p[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
--p[i];
}
if (p[i] < 0) {
// Skip this point
convertedP.clear();
continue;
}
convertedP.push_back(storm::utility::convertNumber<ValueType>(p[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor);
}
}
if (!convertedP.empty()) {
result.push_back(convertedP);
}
}
// 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]) && singleQuantileValues[i].first > 0) {
std::vector<ValueType> newResultPoint(2);
newResultPoint[i] = storm::utility::convertNumber<ValueType>(singleQuantileValues[i].first - 1) * singleQuantileValues[i].second;
newResultPoint[1-i] = storm::utility::infinity<ValueType>();
result.push_back(newResultPoint);
}
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);
}
filterDominatedPoints(result, optimizationDirections);
} else if (infInfProbSatisfiesFormula) {
// then also zeroZeroProb satisfies the formula, i.e., all bound values satisfy the formula
if (storm::solver::minimize(optimizationDirections[0])) {
@ -378,7 +328,129 @@ namespace storm {
}
template<typename ModelType>
std::pair<uint64_t, typename ModelType::ValueType> QuantileHelper<ModelType>::computeQuantileForDimension(Environment const& env, uint64_t dimension) const {
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;
std::cout << "Starting quantile exploration with: (" << currentEpochValues[0] << ", " << currentEpochValues[1] << ")." << std::endl;
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));
exploredEpochs.insert(epoch);
}
}
ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
std::cout << "Numeric result for epoch " << rewardUnfolding.getEpochManager().toString(currEpoch) << " is " << currValue << std::endl;
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()) {
std::cout << "Found another result point: (" << convertedPoint[0] << ", " << convertedPoint[1] << ")." << std::endl;
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
@ -391,7 +463,6 @@ namespace storm {
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));
@ -416,7 +487,18 @@ namespace storm {
}
ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
std::cout << "Numeric result for epoch " << rewardUnfolding.getEpochManager().toString(currEpoch) << " is " << currValue << std::endl;
bool propertySatisfied = transformedFormula->getBound().isSatisfied(currValue);
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) {
@ -438,7 +520,28 @@ namespace storm {
}
++epochValue;
}
return {epochValue, rewardUnfolding.getDimension(dimension).scalingFactor};
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);
}
} else {
return probabilityBound.isSatisfied(numericResult);
}
}
}
template<typename ModelType>
@ -463,7 +566,11 @@ namespace storm {
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();
@ -482,7 +589,7 @@ namespace storm {
ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch);
STORM_LOG_TRACE("Limit probability for infinity dimensions " << infDimensions << " is " << numericResult << ".");
return transformedFormula->getBound().isSatisfied(numericResult);
return numericResult;
}
template class QuantileHelper<storm::models::sparse::Mdp<double>>;

19
src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h

@ -1,6 +1,7 @@
#pragma once
#include "storm/logic/QuantileFormula.h"
#include "boost/optional.hpp"
namespace storm {
class Environment;
@ -22,16 +23,24 @@ namespace storm {
std::vector<std::vector<ValueType>> computeQuantile(Environment const& env);
private:
std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment const& env);
std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env);
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
std::pair<uint64_t, typename ModelType::ValueType> computeQuantileForDimension(Environment const& env, uint64_t dim) 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
Loading…
Cancel
Save