diff --git a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp index 0df580eab..79eae73bb 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp +++ b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp @@ -88,7 +88,7 @@ namespace storm { if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { std::vector cdfEntry; for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) { - uint64_t offset = rewardUnfolding.getDimension(i).isUpperBounded ? 0 : 1; + uint64_t offset = rewardUnfolding.getDimension(i).boundType == helper::rewardbounded::DimensionBoundType::LowerBound ? 1 : 0; cdfEntry.push_back(storm::utility::convertNumber(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor); } auto const& solution = rewardUnfolding.getInitialStateResult(epoch); diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 28e799691..ff7cb5c98 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -129,7 +129,7 @@ namespace storm { if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { std::vector cdfEntry; for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) { - uint64_t offset = rewardUnfolding.getDimension(i).isUpperBounded ? 0 : 1; + uint64_t offset = rewardUnfolding.getDimension(i).boundType == helper::rewardbounded::DimensionBoundType::LowerBound ? 1 : 0; cdfEntry.push_back(storm::utility::convertNumber(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor); } cdfEntry.push_back(rewardUnfolding.getInitialStateResult(epoch)); diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index b28bcd6f0..3aff389fa 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -127,7 +127,7 @@ namespace storm { if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { std::vector cdfEntry; for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) { - uint64_t offset = rewardUnfolding.getDimension(i).isUpperBounded ? 0 : 1; + uint64_t offset = rewardUnfolding.getDimension(i).boundType == helper::rewardbounded::DimensionBoundType::LowerBound ? 1 : 0; cdfEntry.push_back(storm::utility::convertNumber(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor); } cdfEntry.push_back(rewardUnfolding.getInitialStateResult(epoch)); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h index 34a3bef29..6049785ec 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -10,6 +10,14 @@ namespace storm { namespace helper { namespace rewardbounded { + + enum class DimensionBoundType { + Unbounded, // i.e., >=0 or <=B where B approaches infinity + UpperBound, // i.e., <=B where B is either a constant or a variable + LowerBound, // i.e., >B, where B is either a constant or a variable + LowerBoundInfinity // i.e., >B, where B approaches infinity + }; + template struct Dimension { /// The formula describing this dimension @@ -21,11 +29,8 @@ namespace storm { /// A label that indicates the states where this dimension is still relevant (i.e., it is yet unknown whether the corresponding objective holds) boost::optional memoryLabel; - /// True iff the objective is bounded with either an upper bound or a lower bound which is not >= 0 - bool isBounded; - - /// True iff the objective is upper bounded, false if it has a lower bound or no bound at all. - bool isUpperBounded; + /// The type of the bound on this dimension. + DimensionBoundType boundType; /// Multiplying an epoch value with this factor yields the reward/cost in the original domain. ValueType scalingFactor; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 1d57dd2b7..450b04b8c 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -12,6 +12,7 @@ #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/Dtmc.h" +#include "storm/storage/expressions/Expressions.h" #include "storm/transformer/EndComponentEliminator.h" @@ -31,8 +32,9 @@ namespace storm { } template - MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::function&, EpochManager const&)> const& epochStepsCallback) : model(model) { - + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::set const& infinityBoundVariables) : model(model) { + STORM_LOG_TRACE("initializing multi-dimensional reward unfolding for formula " << *objectiveFormula << " and " << infinityBoundVariables.size() << " bound variables should approach infinity."); + if (objectiveFormula->isProbabilityOperatorFormula()) { if (objectiveFormula->getSubformula().isMultiObjectiveFormula()) { for (auto const& subFormula : objectiveFormula->getSubformula().asMultiObjectiveFormula().getSubformulas()) { @@ -53,33 +55,29 @@ namespace storm { objective.considersComplementaryEvent = false; objectives.push_back(std::move(objective)); - initialize(epochStepsCallback); + initialize(infinityBoundVariables); } template - void MultiDimensionalRewardUnfolding::initialize(std::function&, EpochManager const&)> const& epochStepsCallback) { + void MultiDimensionalRewardUnfolding::initialize(std::set const& infinityBoundVariables) { maxSolutionsStored = 0; STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives."); std::vector epochSteps; - initializeObjectives(epochSteps); + initializeObjectives(epochSteps, infinityBoundVariables); + initializeMemoryProduct(epochSteps); - if (epochStepsCallback) { - epochStepsCallback(epochSteps, epochManager); - } // collect which epoch steps are possible possibleEpochSteps.clear(); for (auto const& step : epochSteps) { possibleEpochSteps.insert(step); } - - initializeMemoryProduct(epochSteps); } template - void MultiDimensionalRewardUnfolding::initializeObjectives(std::vector& epochSteps) { + void MultiDimensionalRewardUnfolding::initializeObjectives(std::vector& epochSteps, std::set const& infinityBoundVariables) { std::vector> dimensionWiseEpochSteps; // collect the time-bounded subobjectives for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { @@ -99,28 +97,40 @@ namespace storm { // for simplicity we do not allow interval formulas. STORM_LOG_THROW(!subformula.hasLowerBound(dim) || !subformula.hasUpperBound(dim), storm::exceptions::NotSupportedException, "Bounded until formulas are only supported by this method if they consider either an upper bound or a lower bound. Got " << subformula << " instead."); // lower bounded until formulas with non-trivial left hand side are excluded as this would require some additional effort (in particular the ProductModel::transformMemoryState method). - STORM_LOG_THROW(dimension.isUpperBounded || subformula.getLeftSubformula(dim).isTrueFormula(), storm::exceptions::NotSupportedException, "Lower bounded until formulas are only supported by this method if the left subformula is 'true'. Got " << subformula << " instead."); - dimension.isUpperBounded = subformula.hasUpperBound(dim); + STORM_LOG_THROW(subformula.hasUpperBound(dim) || subformula.getLeftSubformula(dim).isTrueFormula(), storm::exceptions::NotSupportedException, "Lower bounded until formulas are only supported by this method if the left subformula is 'true'. Got " << subformula << " instead."); + // Treat formulas that aren't acutally bounded differently - if ((!subformula.hasLowerBound(dim) && !subformula.hasUpperBound(dim)) || (subformula.hasLowerBound(dim) && !subformula.isLowerBoundStrict(dim) && !subformula.getLowerBound(dim).containsVariables() && storm::utility::isZero(subformula.getLowerBound(dim).evaluateAsRational()))) { + bool formulaUnbounded = (!subformula.hasLowerBound(dim) && !subformula.hasUpperBound(dim)) + || (subformula.hasLowerBound(dim) && !subformula.isLowerBoundStrict(dim) && !subformula.getLowerBound(dim).containsVariables() && storm::utility::isZero(subformula.getLowerBound(dim).evaluateAsRational())) + || (subformula.hasUpperBound(dim) && subformula.getUpperBound(dim).isVariable() && infinityBoundVariables.count(subformula.getUpperBound(dim).getBaseExpression().asVariableExpression().getVariable()) > 0); + if (formulaUnbounded) { dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 0)); dimension.scalingFactor = storm::utility::zero(); - dimension.isBounded = false; - } else if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { - dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); - dimension.scalingFactor = storm::utility::one(); - dimension.isBounded = true; + dimension.boundType = DimensionBoundType::Unbounded; } else { - STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); - std::string const& rewardName = subformula.getTimeBoundReference(dim).getRewardName(); - STORM_LOG_THROW(this->model.hasRewardModel(rewardName), storm::exceptions::IllegalArgumentException, "No reward model with name '" << rewardName << "' found."); - auto const& rewardModel = this->model.getRewardModel(rewardName); - STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Transition rewards are currently not supported as reward bounds."); - std::vector actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix()); - auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(actionRewards); - dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first)); - dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second); - dimension.isBounded = true; + if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { + dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); + dimension.scalingFactor = storm::utility::one(); + } else { + STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); + std::string const& rewardName = subformula.getTimeBoundReference(dim).getRewardName(); + STORM_LOG_THROW(this->model.hasRewardModel(rewardName), storm::exceptions::IllegalArgumentException, "No reward model with name '" << rewardName << "' found."); + auto const& rewardModel = this->model.getRewardModel(rewardName); + STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Transition rewards are currently not supported as reward bounds."); + std::vector actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix()); + auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(actionRewards); + dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first)); + dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second); + } + if (subformula.hasLowerBound(dim)) { + if (subformula.getLowerBound(dim).isVariable() && infinityBoundVariables.count(subformula.getLowerBound(dim).getBaseExpression().asVariableExpression().getVariable()) > 0) { + dimension.boundType = DimensionBoundType::LowerBoundInfinity; + } else { + dimension.boundType = DimensionBoundType::LowerBound; + } + } else { + dimension.boundType = DimensionBoundType::UpperBound; + } } dimensions.emplace_back(std::move(dimension)); } @@ -129,9 +139,9 @@ namespace storm { for (uint64_t dim = 0; dim < subformula.getDimension(); ++dim) { Dimension dimension; dimension.formula = subformula.restrictToDimension(dim); + STORM_LOG_THROW(!(dimension.formula->asCumulativeRewardFormula().getBound().isVariable() && infinityBoundVariables.count(dimension.formula->asCumulativeRewardFormula().getBound().getBaseExpression().asVariableExpression().getVariable()) > 0), storm::exceptions::NotSupportedException, "Letting cumulative reward bounds approach infinite is not supported."); dimension.objectiveIndex = objIndex; - dimension.isUpperBounded = true; - dimension.isBounded = true; + dimension.boundType = DimensionBoundType::UpperBound; if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); dimension.scalingFactor = storm::utility::one(); @@ -168,7 +178,7 @@ namespace storm { objDimensions.set(currDim); } for (uint64_t currDim = dim; currDim < dim + objDimensionCount; ++currDim ) { - if (!objDimensionsCanBeSatisfiedIndividually || dimensions[currDim].isUpperBounded) { + if (!objDimensionsCanBeSatisfiedIndividually || dimensions[currDim].boundType == DimensionBoundType::UpperBound) { dimensions[currDim].dependentDimensions = objDimensions; } else { dimensions[currDim].dependentDimensions = storm::storage::BitVector(dimensions.size(), false); @@ -198,7 +208,7 @@ namespace storm { // Set the maximal values we need to consider for each dimension computeMaxDimensionValues(); - + translateLowerBoundInfinityDimensions(epochSteps); } template @@ -213,7 +223,7 @@ namespace storm { bool isStrict = false; storm::logic::Formula const& dimFormula = *dimensions[dim].formula; if (dimFormula.isBoundedUntilFormula()) { - assert(!dimFormula.asBoundedUntilFormula().isMultiDimensional()); + assert(!dimFormula.asBoundedUntilFormula().hasMultiDimensionalSubformulas()); if (dimFormula.asBoundedUntilFormula().hasUpperBound()) { STORM_LOG_ASSERT(!dimFormula.asBoundedUntilFormula().hasLowerBound(), "Bounded until formulas with interval bounds are not supported."); bound = dimFormula.asBoundedUntilFormula().getUpperBound(); @@ -232,11 +242,11 @@ namespace storm { if (!bound.containsVariables()) { // We always consider upper bounds to be non-strict and lower bounds to be strict. // Thus, >=N would become >N-1. However, note that the case N=0 is treated separately. - if (dimensions[dim].isBounded) { + if (dimensions[dim].boundType == DimensionBoundType::LowerBound || dimensions[dim].boundType == DimensionBoundType::UpperBound) { ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); discretizedBound /= dimensions[dim].scalingFactor; if (storm::utility::isInteger(discretizedBound)) { - if (isStrict == dimensions[dim].isUpperBounded) { + if (isStrict == (dimensions[dim].boundType == DimensionBoundType::UpperBound)) { discretizedBound -= storm::utility::one(); } } else { @@ -250,14 +260,68 @@ namespace storm { } } + template + void MultiDimensionalRewardUnfolding::translateLowerBoundInfinityDimensions(std::vector& epochSteps) { + // Translate lowerBoundedByInfinity dimensions to finite bounds + storm::storage::BitVector infLowerBoundedDimensions(dimensions.size(), false); + storm::storage::BitVector upperBoundedDimensions(dimensions.size(), false); + for (uint64_t dim = 0; dim < dimensions.size(); ++dim) { + infLowerBoundedDimensions.set(dim, dimensions[dim].boundType == DimensionBoundType::LowerBoundInfinity); + upperBoundedDimensions.set(dim, dimensions[dim].boundType == DimensionBoundType::UpperBound); + } + if (!infLowerBoundedDimensions.empty()) { + // We can currently only handle this case for single maximizing bounded until probability objectives. + // The approach is to erase all epochSteps that are not part of an end component and then change the reward bound to '> 0'. + // Then, reaching a reward means reaching an end component where arbitrarily many reward can be collected. + 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(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::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. + for (uint64_t choiceIndex = 0; choiceIndex < model.getNumberOfChoices(); ++choiceIndex) { + for (auto const& dim : upperBoundedDimensions) { + if (epochManager.getDimensionOfEpoch(epochSteps[choiceIndex], dim) != 0) { + choicesWithoutUpperBoundedStep.set(choiceIndex, false); + break; + } + } + } + } + storm::storage::MaximalEndComponentDecomposition mecDecomposition(model.getTransitionMatrix(), model.getBackwardTransitions(), storm::storage::BitVector(model.getNumberOfStates(), true), choicesWithoutUpperBoundedStep); + storm::storage::BitVector nonMecChoices(model.getNumberOfChoices(), true); + for (auto const& mec : mecDecomposition) { + for (auto const& stateChoicesPair : mec) { + for (auto const& choice : stateChoicesPair.second) { + nonMecChoices.set(choice, false); + } + } + } + for (auto const& choice : nonMecChoices) { + for (auto const& dim : infLowerBoundedDimensions) { + epochManager.setDimensionOfEpoch(epochSteps[choice], dim, 0); + } + } + + // Translate the dimension to '>0' + for (auto const& dim : infLowerBoundedDimensions) { + dimensions[dim].boundType = DimensionBoundType::LowerBound; + dimensions[dim].maxValue = 0; + } + } + } + template typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch(bool setUnknownDimsToBottom) { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (dimensions[dim].isBounded && dimensions[dim].maxValue) { + if (dimensions[dim].maxValue) { epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); } else { - STORM_LOG_THROW(setUnknownDimsToBottom || !dimensions[dim].isBounded, storm::exceptions::UnexpectedException, "Tried to obtain the start epoch although not all dimensions are known."); + STORM_LOG_THROW(setUnknownDimsToBottom || dimensions[dim].boundType == DimensionBoundType::Unbounded, storm::exceptions::UnexpectedException, "Tried to obtain the start epoch although no bound on dimension " << dim << " is known."); epochManager.setBottomDimension(startEpoch, dim); } } @@ -318,7 +382,7 @@ namespace storm { bool containsLowerBoundedObjective = false; for (auto const& dimension : dimensions) { - if (!dimension.isUpperBounded) { + if (dimension.boundType == DimensionBoundType::LowerBound) { containsLowerBoundedObjective = true; break; } @@ -348,7 +412,7 @@ namespace storm { bool rewardEarned = !storm::utility::isZero(epochModel.objectiveRewards[objIndex][reducedChoice]); if (rewardEarned) { for (auto const& dim : objectiveDimensions[objIndex]) { - if (dimensions[dim].isUpperBounded == epochManager.isBottomDimension(successorEpoch, dim) && productModel->getMemoryStateManager().isRelevantDimension(memoryState, dim)) { + if ((dimensions[dim].boundType == DimensionBoundType::UpperBound) == epochManager.isBottomDimension(successorEpoch, dim) && productModel->getMemoryStateManager().isRelevantDimension(memoryState, dim)) { rewardEarned = false; break; } @@ -430,7 +494,7 @@ namespace storm { // redirect transitions for the case where the lower reward bounds are not met yet storm::storage::BitVector violatedLowerBoundedDimensions(dimensions.size(), false); for (uint64_t dim = 0; dim < dimensions.size(); ++dim) { - if (!dimensions[dim].isUpperBounded && !epochManager.isBottomDimensionEpochClass(epochClass, dim)) { + if (dimensions[dim].boundType == DimensionBoundType::LowerBound && !epochManager.isBottomDimensionEpochClass(epochClass, dim)) { violatedLowerBoundedDimensions.set(dim); } } @@ -569,6 +633,18 @@ namespace storm { } + template + template::type> + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::getOneSolution() const { + return storm::utility::one; + } + + template + template::type> + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::getOneSolution() const { + return SolutionType(objectives.size(), storm::utility::one()); + } + template template::type> typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const { @@ -775,11 +851,9 @@ namespace storm { typename MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::getStateSolution(Epoch const& epoch, uint64_t const& productState) { auto epochSolutionIt = epochSolutions.find(epoch); STORM_LOG_ASSERT(epochSolutionIt != epochSolutions.end(), "Requested unexisting solution for epoch " << epochManager.toString(epoch) << "."); - auto const& epochSolution = epochSolutionIt->second; - STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution for epoch " << epochManager.toString(epoch) << " at an unexisting product state."); - STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch " << epochManager.toString(epoch) << " at a state for which no solution was stored."); - return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]]; + return getStateSolution(epochSolutionIt->second, productState); } + template typename MultiDimensionalRewardUnfolding::EpochSolution const& MultiDimensionalRewardUnfolding::getEpochSolution(std::map const& solutions, Epoch const& epoch) { auto epochSolutionIt = solutions.find(epoch); @@ -791,6 +865,10 @@ namespace storm { 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."); + 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(); + } return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]]; } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h index 77f1ec64b..6513b1939 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -36,7 +36,15 @@ namespace storm { * */ MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::vector> const& objectives); - MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::function&, EpochManager const&)> const& epochStepsCallback = nullptr); + + /*! + * Initializes the reward unfolding with just a single objective. + * + * @param model The input model + * @param objectiveFormula the formula + * @param infinityBoundVariables if non-empty, reward bounds with these variables are assumed to approach infinity + */ + MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::set const& infinityBoundVariables = {}); ~MultiDimensionalRewardUnfolding() = default; @@ -73,13 +81,20 @@ namespace storm { private: void setCurrentEpochClass(Epoch const& epoch); - void initialize(std::function&, EpochManager const&)> const& epochStepsCallback = nullptr); + void initialize(std::set const& infinityBoundVariables = {}); - void initializeObjectives(std::vector& epochSteps); + void initializeObjectives(std::vector& epochSteps, std::set const& infinityBoundVariables); void computeMaxDimensionValues(); + void translateLowerBoundInfinityDimensions(std::vector& epochSteps); void initializeMemoryProduct(std::vector const& epochSteps); + // Returns a solution only consisting of one + template::type = 0> + SolutionType getOneSolution() const; + template::type = 0> + SolutionType getOneSolution() const; + template::type = 0> SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) 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 57fbf4980..4ed49fa1a 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()) { + 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) { for (uint64_t dim = 0; dim < dimensions.size(); ++dim) { if (!dimensions[dim].memoryLabel) { @@ -112,7 +112,7 @@ namespace storm { bool objectiveContainsLowerBound = false; for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) { - if (!dimensions[globalDimensionIndex].isUpperBounded) { + if (dimensions[globalDimensionIndex].boundType == DimensionBoundType::LowerBound) { objectiveContainsLowerBound = true; break; } @@ -168,8 +168,12 @@ 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 is not supported as we can not reduce this (easily) to an expected reward computation. - STORM_LOG_THROW(!memStatePrimeBV.empty() || initialTransitionStates.empty() || objectiveContainsLowerBound || initialTransitionStates.isDisjointFrom(constraintStates), storm::exceptions::NotSupportedException, "The objective " << *objectives[objIndex].formula << " is already satisfied in an initial state. This special case is not supported."); + // 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); + } + for (auto const& initState : initialTransitionStates) { objMemoryBuilder.setInitialMemoryState(initState, memStatePrime); } @@ -231,14 +235,34 @@ namespace storm { } // Initialize the reachable states with the initial states - EpochClass initEpochClass = epochManager.getEpochClass(epochManager.getZeroEpoch()); - auto memStateIt = memory.getInitialMemoryStates().begin(); - for (auto const& initState : model.getInitialStates()) { - uint64_t transformedMemoryState = transformMemoryState(memoryStateMap[*memStateIt], initEpochClass, memoryStateManager.getInitialMemoryState()); - reachableProductStates[transformedMemoryState].set(initState, true); - ++memStateIt; + // If the bound is not known (e.g., when computing quantiles) we don't know the initial epoch class. Hence, we consider all possibilities. + std::vector initEpochClasses; + initEpochClasses.push_back(epochManager.getEpochClass(epochManager.getZeroEpoch())); + for (uint64_t dim = 0; dim < dimensions.size(); ++dim) { + Dimension const& dimension = dimensions[dim]; + if (dimension.boundType == DimensionBoundType::Unbounded) { + // For unbounded dimensions we are only interested in the bottom class. + for (auto& ec : initEpochClasses) { + epochManager.setDimensionOfEpochClass(ec, dim, true); + } + } else if (!dimension.maxValue) { + // If no max value is known, we have to consider all possibilities + std::vector newEcs = initEpochClasses; + for (auto& ec : newEcs) { + epochManager.setDimensionOfEpochClass(ec, dim, true); + } + initEpochClasses.insert(initEpochClasses.end(), newEcs.begin(), newEcs.end()); + } + } + for (auto const& initEpochClass : initEpochClasses) { + auto memStateIt = memory.getInitialMemoryStates().begin(); + for (auto const& initState : model.getInitialStates()) { + uint64_t transformedMemoryState = transformMemoryState(memoryStateMap[*memStateIt], initEpochClass, memoryStateManager.getInitialMemoryState()); + reachableProductStates[transformedMemoryState].set(initState, true); + ++memStateIt; + } + assert(memStateIt == memory.getInitialMemoryStates().end()); } - assert(memStateIt == memory.getInitialMemoryStates().end()); // Find the reachable epoch classes std::set possibleSteps(originalModelSteps.begin(), originalModelSteps.end()); @@ -303,7 +327,7 @@ namespace storm { template uint64_t ProductModel::getProductState(uint64_t const& modelState, MemoryState const& memoryState) const { - STORM_LOG_ASSERT(productStateExists(modelState, memoryState), "Tried to obtain a state in the model-memory-product that does not exist"); + STORM_LOG_ASSERT(productStateExists(modelState, memoryState), "Tried to obtain state (" << modelState << ", " << memoryStateManager.toString(memoryState) << ") in the model-memory-product which does not exist"); return modelMemoryToProductStateMap[modelState * memoryStateManager.getUpperMemoryStateBound() + memoryState]; } @@ -369,7 +393,7 @@ namespace storm { // find out whether objective reward should be earned within this epoch class bool collectRewardInEpoch = true; for (auto const& subObjIndex : relevantObjectives) { - if (dimensions[dimensionIndexMap[subObjIndex]].isUpperBounded && epochManager.isBottomDimensionEpochClass(epochClass, dimensionIndexMap[subObjIndex])) { + if (dimensions[dimensionIndexMap[subObjIndex]].boundType == DimensionBoundType::UpperBound && epochManager.isBottomDimensionEpochClass(epochClass, dimensionIndexMap[subObjIndex])) { collectRewardInEpoch = false; break; } @@ -487,9 +511,9 @@ namespace storm { } } - // Also treat dimensions without a priori bound + // Also treat dimensions without a priori bound. Unbounded dimensions need no further treatment as for these only the 'bottom' class is relevant. for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (!dimensions[dim].maxValue) { + if (dimensions[dim].boundType != DimensionBoundType::Unbounded && !dimensions[dim].maxValue) { std::vector newClasses; for (auto const& c : reachableEpochClasses) { auto newClass = c; @@ -511,7 +535,7 @@ namespace storm { for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { if (epochManager.isBottomDimensionEpochClass(epochClass, dim)) { bottomDimensions.set(dim, true); - if (dimensions[dim].isBounded && dimensions[dim].maxValue) { + if (dimensions[dim].boundType != DimensionBoundType::Unbounded && dimensions[dim].maxValue) { considerInitialStates = false; } } @@ -602,7 +626,7 @@ namespace storm { for (auto const& dim : objDimensions) { auto const& dimension = dimensions[dim]; if (dimension.memoryLabel) { - bool dimUpperBounded = dimension.isUpperBounded; + bool dimUpperBounded = dimension.boundType == DimensionBoundType::UpperBound; bool dimBottom = epochManager.isBottomDimensionEpochClass(epochClass, dim); if (dimUpperBounded && dimBottom && memoryStateManager.isRelevantDimension(predecessorMemoryState, dim)) { STORM_LOG_ASSERT(objDimensions == dimension.dependentDimensions, "Unexpected set of dependent dimensions"); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h index 3f5b2a9a7..92036761e 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h @@ -46,7 +46,9 @@ 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(); private: @@ -75,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 }; } } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index a06907e2f..84a28a68c 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -35,7 +35,9 @@ 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 - + // 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?) @@ -56,18 +58,20 @@ namespace storm { std::vector> leftSubformulas, rightSubformulas; std::vector> lowerBounds, upperBounds; std::vector timeBoundReferences; - + for (uint64_t dim = 0; dim < origBoundedUntil.getDimension(); ++dim) { - leftSubformulas.push_back(origBoundedUntil.getLeftSubformula(dim).asSharedPointer()); - rightSubformulas.push_back(origBoundedUntil.getRightSubformula(dim).asSharedPointer()); + if (origBoundedUntil.hasMultiDimensionalSubformulas()) { + leftSubformulas.push_back(origBoundedUntil.getLeftSubformula(dim).asSharedPointer()); + rightSubformulas.push_back(origBoundedUntil.getRightSubformula(dim).asSharedPointer()); + } timeBoundReferences.push_back(origBoundedUntil.getTimeBoundReference(dim)); if (transformations[dim] == BoundTransformation::None) { - if (origBoundedUntil.hasLowerBound()) { + if (origBoundedUntil.hasLowerBound(dim)) { lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); } else { lowerBounds.push_back(boost::none); } - if (origBoundedUntil.hasUpperBound()) { + if (origBoundedUntil.hasUpperBound(dim)) { upperBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim))); } else { upperBounds.push_back(boost::none); @@ -91,7 +95,12 @@ namespace storm { } } } - auto newBoundedUntil = std::make_shared(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences); + std::shared_ptr newBoundedUntil; + if (origBoundedUntil.hasMultiDimensionalSubformulas()) { + newBoundedUntil = std::make_shared(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences); + } else { + newBoundedUntil = std::make_shared(origBoundedUntil.getLeftSubformula().asSharedPointer(), origBoundedUntil.getRightSubformula().asSharedPointer(), lowerBounds, upperBounds, timeBoundReferences); + } return std::make_shared(newBoundedUntil, boundedUntilOperator.getOperatorInformation()); } @@ -103,17 +112,19 @@ namespace storm { template storm::storage::BitVector QuantileHelper::getOpenDimensions() const { auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); - storm::storage::BitVector res(getDimension()); + storm::storage::BitVector res(getDimension(), false); for (uint64_t dim = 0; dim < getDimension(); ++dim) { - res.set(dim, boundedUntil.hasLowerBound(dim) ? boundedUntil.getLowerBound(dim).containsVariables() : boundedUntil.getUpperBound(dim).containsVariables()); + auto const& bound = boundedUntil.hasLowerBound(dim) ? boundedUntil.getLowerBound(dim) : boundedUntil.getUpperBound(dim); + if (bound.containsVariables()) { + res.set(dim, true); + } } return res; } template storm::solver::OptimizationDirection const& QuantileHelper::getOptimizationDirForDimension(uint64_t const& dim) const { - auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); - storm::expressions::Variable const& dimVar = (boundedUntil.hasLowerBound() ? boundedUntil.getLowerBound(dim) : boundedUntil.getUpperBound(dim)).getBaseExpression().asVariableExpression().getVariable(); + storm::expressions::Variable const& dimVar = getVariableForDimension(dim); for (auto const& boundVar : quantileFormula.getBoundVariables()) { if (boundVar.second == dimVar) { return boundVar.first; @@ -122,6 +133,12 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "The bound variable '" << dimVar.getName() << "' is not specified within the quantile formula '" << quantileFormula << "'."); return quantileFormula.getOptimizationDirection(); } + + template + storm::expressions::Variable const& QuantileHelper::getVariableForDimension(uint64_t const& dim) const { + auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); + return (boundedUntil.hasLowerBound(dim) ? boundedUntil.getLowerBound(dim) : boundedUntil.getUpperBound(dim)).getBaseExpression().asVariableExpression().getVariable(); + } template storm::storage::BitVector QuantileHelper::getDimensionsForVariable(storm::expressions::Variable const& var) const { @@ -145,13 +162,13 @@ namespace storm { uint64_t dimension = *getOpenDimensions().begin(); auto const& boundedUntilOperator = quantileFormula.getSubformula().asProbabilityOperatorFormula(); - bool maxProbSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeExtremalValue(env, storm::storage::BitVector(getDimension(), false))); - bool minProbSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeExtremalValue(env, getOpenDimensions())); - if (maxProbSatisfiesFormula != minProbSatisfiesFormula) { + bool zeroSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeLimitValue(env, storm::storage::BitVector(getDimension(), false))); + bool infSatisfiesFormula = boundedUntilOperator.getBound().isSatisfied(computeLimitValue(env, getOpenDimensions())); + if (zeroSatisfiesFormula != infSatisfiesFormula) { auto quantileRes = computeQuantileForDimension(env, dimension); result = {{storm::utility::convertNumber(quantileRes.first) * quantileRes.second}}; - } else if (maxProbSatisfiesFormula) { - // i.e., all bound values satisfy the formula + } 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() }}; } else { @@ -223,31 +240,26 @@ 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 maxmaxProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, storm::storage::BitVector(getDimension(), false))); - bool minminProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1])); - if (maxmaxProbSatisfiesFormula != minminProbSatisfiesFormula) { + bool infInfProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeLimitValue(env, storm::storage::BitVector(getDimension(), false))); + bool zeroZeroProbSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeLimitValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1])); + if (infInfProbSatisfiesFormula != zeroZeroProbSatisfiesFormula) { std::vector> singleQuantileValues; std::vector> resultPoints; const int64_t infinity = std::numeric_limits().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 - uint64_t indexToMinimizeProb = lowerCostBounds[i] ? (1-i) : i; - bool zeroInfSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeExtremalValue(env, dimensionsAsBitVector[indexToMinimizeProb])); + bool zeroInfSatisfiesFormula = probOpFormula.getBound().isSatisfied(computeLimitValue(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 singleQuantileValues.emplace_back(0, storm::utility::zero()); } else { - // Compute quantile where 1-i approaches infinity + // 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])); std::cout << ".. Result is " << singleQuantileValues.back().first << std::endl; - // When maximizing bounds, the computed quantile value is sat for all values of the other bound. if (!storm::solver::minimize(optimizationDirections[i])) { - std::vector newResultPoint(2); - newResultPoint[i] = singleQuantileValues.back().first; - newResultPoint[1-i] = infinity; - resultPoints.push_back(newResultPoint); + // 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; } @@ -334,9 +346,19 @@ namespace storm { 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 newResultPoint(2); + newResultPoint[i] = storm::utility::convertNumber(singleQuantileValues[i].first - 1) * singleQuantileValues[i].second; + newResultPoint[1-i] = storm::utility::infinity(); + result.push_back(newResultPoint); + } + } filterDominatedPoints(result, optimizationDirections); - } else if (maxmaxProbSatisfiesFormula) { - // i.e., all bound values satisfy the formula + + } 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 { @@ -360,49 +382,19 @@ namespace storm { // 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::vector bts(getDimension(), BoundTransformation::None); - storm::storage::BitVector otherLowerBoundedDimensions(getDimension(), false); + std::set infinityVariables; for (auto const& d : getOpenDimensions()) { if (d != dimension) { - if (quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasLowerBound(d)) { - bts[d] = BoundTransformation::GreaterZero; - otherLowerBoundedDimensions.set(d, true); - } else { - bts[d] = BoundTransformation::GreaterEqualZero; - } + 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)); - storm::storage::BitVector nonMecChoices; - if (!otherLowerBoundedDimensions.empty()) { - STORM_LOG_ASSERT(!upperCostBound, "It is not possible to let other open dimensions approach infinity if this is an upper cost bound and others are lower cost bounds."); - // Get the choices that do not lie on a mec - nonMecChoices.resize(model.getNumberOfChoices(), true); - auto mecDecomposition = storm::storage::MaximalEndComponentDecomposition(model); - for (auto const& mec : mecDecomposition) { - for (auto const& stateChoicesPair : mec) { - for (auto const& choice : stateChoicesPair.second) { - nonMecChoices.set(choice, false); - } - } - } - } - - auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts); - MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, [&](std::vector& epochSteps, EpochManager const& epochManager) { - if (!otherLowerBoundedDimensions.empty()) { - for (auto const& choice : nonMecChoices) { - for (auto const& dim : otherLowerBoundedDimensions) { - epochManager.setDimensionOfEpoch(epochSteps[choice], dim, 0); - } - } - } - }); - - // initialize data that will be needed for each epoch auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); auto upperBound = rewardUnfolding.getUpperObjectiveBound(); @@ -450,58 +442,27 @@ namespace storm { } template - typename ModelType::ValueType QuantileHelper::computeExtremalValue(Environment const& env, storm::storage::BitVector const& minimizingDimensions) const { - // For maximizing in a dimension, we can simply 'drop' the bound by replacing it with >=0 - // For minimizing an upper-bounded dimension, we can replace it with <=0 - // For minimizing a lower-bounded dimension, the lower bound needs to approach infinity. - // To compute this limit probability, we only consider epoch steps that lie in an end component and check for the bound >0 instead. - // Notice, however, that this approach fails if we try to minimize for a lower and an upper bounded dimension - + 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); - storm::storage::BitVector minimizingLowerBoundedDimensions(getDimension(), false); - + std::set infinityVariables; for (auto const& d : getOpenDimensions()) { - if (minimizingDimensions.get(d)) { - bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d); + bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d); + if (infDimensions.get(d)) { + infinityVariables.insert(getVariableForDimension(d)); + } else { if (upperCostBound) { bts[d] = BoundTransformation::LessEqualZero; - STORM_LOG_ASSERT(std::find(bts.begin(), bts.end(), BoundTransformation::GreaterZero) == bts.end(), "Unable to compute extremal value for minimizing lower and upper bounded dimensions."); } else { - bts[d] = BoundTransformation::GreaterZero; - minimizingLowerBoundedDimensions.set(d, true); - STORM_LOG_ASSERT(std::find(bts.begin(), bts.end(), BoundTransformation::LessEqualZero) == bts.end(), "Unable to compute extremal value for minimizing lower and upper bounded dimensions."); + bts[d] = BoundTransformation::GreaterEqualZero; } - } else { - bts[d] = BoundTransformation::GreaterEqualZero; } } - auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts); - storm::storage::BitVector nonMecChoices; - if (!minimizingLowerBoundedDimensions.empty()) { - // Get the choices that do not lie on a mec - nonMecChoices.resize(model.getNumberOfChoices(), true); - auto mecDecomposition = storm::storage::MaximalEndComponentDecomposition(model); - for (auto const& mec : mecDecomposition) { - for (auto const& stateChoicesPair : mec) { - for (auto const& choice : stateChoicesPair.second) { - nonMecChoices.set(choice, false); - } - } - } - } - std::cout << "nonmec choices are " << nonMecChoices << std::endl; - - MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, [&](std::vector& epochSteps, EpochManager const& epochManager) { - if (!minimizingLowerBoundedDimensions.empty()) { - for (auto const& choice : nonMecChoices) { - for (auto const& dim : minimizingLowerBoundedDimensions) { - epochManager.setDimensionOfEpoch(epochSteps[choice], dim, 0); - } - } - } - }); + MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, infinityVariables); // Get lower and upper bounds for the solution. auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); @@ -520,7 +481,7 @@ namespace storm { } ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch); - STORM_LOG_TRACE("Extremal probability for minimizing dimensions " << minimizingDimensions << " is " << numericResult << "."); + STORM_LOG_TRACE("Limit probability for infinity dimensions " << infDimensions << " is " << numericResult << "."); return transformedFormula->getBound().isSatisfied(numericResult); } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h index a0cd22fdb..ce05b5f5a 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -25,10 +25,9 @@ namespace storm { std::vector> computeTwoDimensionalQuantile(Environment const& env); /*! - * Computes the infimum over bound values in minimizing dimensions / the supremum over bound values in maximizing dimensions - * @param minimizingDimensions marks dimensions which should be minimized. The remaining dimensions are either not open or maximizing. + * Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero. */ - ValueType computeExtremalValue(Environment const& env, storm::storage::BitVector const& minimizingDimensions) const; + ValueType computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const; /// Computes the quantile with respect to the given dimension @@ -47,6 +46,7 @@ namespace storm { storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var) const; storm::solver::OptimizationDirection const& getOptimizationDirForDimension(uint64_t const& dim) const; + storm::expressions::Variable const& getVariableForDimension(uint64_t const& dim) const; ModelType const& model; storm::logic::QuantileFormula const& quantileFormula;