From 7e66787c9c8a956629ea934a3e090bd9062a19de Mon Sep 17 00:00:00 2001 From: TimQu Date: Thu, 24 Jan 2019 12:42:57 +0100 Subject: [PATCH 01/20] logic: Added QuantileFormulas. --- src/storm/logic/CloneVisitor.cpp | 5 + src/storm/logic/CloneVisitor.h | 1 + src/storm/logic/Formula.cpp | 4 + src/storm/logic/Formula.h | 1 + src/storm/logic/FormulaInformationVisitor.cpp | 4 + src/storm/logic/FormulaInformationVisitor.h | 1 + src/storm/logic/FormulaVisitor.h | 1 + src/storm/logic/Formulas.h | 1 + src/storm/logic/FormulasForwardDeclarations.h | 1 + src/storm/logic/FragmentChecker.cpp | 11 +++ src/storm/logic/FragmentChecker.h | 1 + src/storm/logic/FragmentSpecification.cpp | 43 ++++++++- src/storm/logic/FragmentSpecification.h | 11 +++ .../LiftableTransitionRewardsVisitor.cpp | 4 + .../logic/LiftableTransitionRewardsVisitor.h | 1 + src/storm/logic/QuantileFormula.cpp | 96 +++++++++++++++++++ src/storm/logic/QuantileFormula.h | 41 ++++++++ src/storm/logic/ToExpressionVisitor.cpp | 4 + src/storm/logic/ToExpressionVisitor.h | 1 + src/storm/storage/jani/JSONExporter.cpp | 4 + src/storm/storage/jani/JSONExporter.h | 1 + 21 files changed, 236 insertions(+), 1 deletion(-) create mode 100644 src/storm/logic/QuantileFormula.cpp create mode 100644 src/storm/logic/QuantileFormula.h diff --git a/src/storm/logic/CloneVisitor.cpp b/src/storm/logic/CloneVisitor.cpp index 4b427b1dc..cb24e7d8b 100644 --- a/src/storm/logic/CloneVisitor.cpp +++ b/src/storm/logic/CloneVisitor.cpp @@ -108,6 +108,11 @@ namespace storm { return std::static_pointer_cast(std::make_shared(subformulas)); } + boost::any CloneVisitor::visit(QuantileFormula const& f, boost::any const& data) const { + std::shared_ptr subformula = boost::any_cast>(f.getSubformula().accept(*this, data)); + return std::static_pointer_cast(std::make_shared(f.getBoundVariables(), subformula)); + } + boost::any CloneVisitor::visit(NextFormula const& f, boost::any const& data) const { std::shared_ptr subformula = boost::any_cast>(f.getSubformula().accept(*this, data)); return std::static_pointer_cast(std::make_shared(subformula)); diff --git a/src/storm/logic/CloneVisitor.h b/src/storm/logic/CloneVisitor.h index c03f2e4fc..4b1d620ee 100644 --- a/src/storm/logic/CloneVisitor.h +++ b/src/storm/logic/CloneVisitor.h @@ -26,6 +26,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const override; virtual boost::any visit(NextFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/Formula.cpp b/src/storm/logic/Formula.cpp index 9dd2b0fba..fe3a1e41a 100644 --- a/src/storm/logic/Formula.cpp +++ b/src/storm/logic/Formula.cpp @@ -23,6 +23,10 @@ namespace storm { return false; } + bool Formula::isQuantileFormula() const { + return false; + } + bool Formula::isBinaryStateFormula() const { return false; } diff --git a/src/storm/logic/Formula.h b/src/storm/logic/Formula.h index 91ccdf342..bd88374a3 100644 --- a/src/storm/logic/Formula.h +++ b/src/storm/logic/Formula.h @@ -48,6 +48,7 @@ namespace storm { virtual bool isUnaryBooleanStateFormula() const; virtual bool isMultiObjectiveFormula() const; + virtual bool isQuantileFormula() const; // Operator formulas. virtual bool isOperatorFormula() const; diff --git a/src/storm/logic/FormulaInformationVisitor.cpp b/src/storm/logic/FormulaInformationVisitor.cpp index 0aab5e04e..b3ad356db 100644 --- a/src/storm/logic/FormulaInformationVisitor.cpp +++ b/src/storm/logic/FormulaInformationVisitor.cpp @@ -93,6 +93,10 @@ namespace storm { return result; } + boost::any FormulaInformationVisitor::visit(QuantileFormula const& f, boost::any const& data) const { + return f.getSubformula().accept(*this, data); + } + boost::any FormulaInformationVisitor::visit(NextFormula const& f, boost::any const& data) const { return boost::any_cast(f.getSubformula().accept(*this, data)).setContainsNextFormula(); } diff --git a/src/storm/logic/FormulaInformationVisitor.h b/src/storm/logic/FormulaInformationVisitor.h index 134ad0bc4..7ba8b70b9 100644 --- a/src/storm/logic/FormulaInformationVisitor.h +++ b/src/storm/logic/FormulaInformationVisitor.h @@ -25,6 +25,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const override; virtual boost::any visit(NextFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/FormulaVisitor.h b/src/storm/logic/FormulaVisitor.h index f0abd0436..1eaf6791d 100644 --- a/src/storm/logic/FormulaVisitor.h +++ b/src/storm/logic/FormulaVisitor.h @@ -26,6 +26,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const = 0; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(NextFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const = 0; diff --git a/src/storm/logic/Formulas.h b/src/storm/logic/Formulas.h index bba8afea0..bdda5c863 100644 --- a/src/storm/logic/Formulas.h +++ b/src/storm/logic/Formulas.h @@ -19,6 +19,7 @@ #include "storm/logic/StateFormula.h" #include "storm/logic/LongRunAverageOperatorFormula.h" #include "storm/logic/MultiObjectiveFormula.h" +#include "storm/logic/QuantileFormula.h" #include "storm/logic/TimeOperatorFormula.h" #include "storm/logic/TotalRewardFormula.h" #include "storm/logic/UnaryBooleanStateFormula.h" diff --git a/src/storm/logic/FormulasForwardDeclarations.h b/src/storm/logic/FormulasForwardDeclarations.h index 3f3f1dfaf..170316fc7 100644 --- a/src/storm/logic/FormulasForwardDeclarations.h +++ b/src/storm/logic/FormulasForwardDeclarations.h @@ -21,6 +21,7 @@ namespace storm { class LongRunAverageOperatorFormula; class LongRunAverageRewardFormula; class MultiObjectiveFormula; + class QuantileFormula; class NextFormula; class OperatorFormula; struct OperatorInformation; diff --git a/src/storm/logic/FragmentChecker.cpp b/src/storm/logic/FragmentChecker.cpp index adc96728b..4bda466c4 100644 --- a/src/storm/logic/FragmentChecker.cpp +++ b/src/storm/logic/FragmentChecker.cpp @@ -27,6 +27,9 @@ namespace storm { if (specification.isMultiObjectiveFormulaAtTopLevelRequired()) { result &= f.isMultiObjectiveFormula(); } + if (specification.isQuantileFormulaAtTopLevelRequired()) { + result &= f.isQuantileFormula(); + } return result; } @@ -229,6 +232,14 @@ namespace storm { return result; } + boost::any FragmentChecker::visit(QuantileFormula const& f, boost::any const& data) const { + InheritedInformation const& inherited = boost::any_cast(data); + if (!inherited.getSpecification().areQuantileFormulasAllowed()) { + return false; + } + return f.getSubformula().accept(*this, data); + } + boost::any FragmentChecker::visit(NextFormula const& f, boost::any const& data) const { InheritedInformation const& inherited = boost::any_cast(data); bool result = inherited.getSpecification().areNextFormulasAllowed(); diff --git a/src/storm/logic/FragmentChecker.h b/src/storm/logic/FragmentChecker.h index cb3ddb95d..485a205aa 100644 --- a/src/storm/logic/FragmentChecker.h +++ b/src/storm/logic/FragmentChecker.h @@ -26,6 +26,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const override; virtual boost::any visit(NextFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/FragmentSpecification.cpp b/src/storm/logic/FragmentSpecification.cpp index 9b5b8b17e..46d7c4075 100644 --- a/src/storm/logic/FragmentSpecification.cpp +++ b/src/storm/logic/FragmentSpecification.cpp @@ -110,6 +110,27 @@ namespace storm { return multiObjective; } + FragmentSpecification quantiles() { + FragmentSpecification quantiles = propositional(); + + quantiles.setQuantileFormulasAllowed(true); + quantiles.setQuantileFormulaAtTopLevelRequired(true); + quantiles.setProbabilityOperatorsAllowed(true); + quantiles.setRewardOperatorsAllowed(true); + quantiles.setBoundedUntilFormulasAllowed(true); + quantiles.setStepBoundedUntilFormulasAllowed(true); + quantiles.setTimeBoundedUntilFormulasAllowed(true); + quantiles.setRewardBoundedUntilFormulasAllowed(true); + quantiles.setStepBoundedCumulativeRewardFormulasAllowed(true); + quantiles.setTimeBoundedCumulativeRewardFormulasAllowed(true); + quantiles.setRewardBoundedCumulativeRewardFormulasAllowed(true); + quantiles.setCumulativeRewardFormulasAllowed(true); + quantiles.setMultiDimensionalBoundedUntilFormulasAllowed(true); + quantiles.setMultiDimensionalCumulativeRewardFormulasAllowed(true); + + return quantiles; + } + FragmentSpecification::FragmentSpecification() { probabilityOperator = false; rewardOperator = false; @@ -117,6 +138,7 @@ namespace storm { longRunAverageOperator = false; multiObjectiveFormula = false; + quantileFormula = false; globallyFormula = false; reachabilityProbabilityFormula = false; @@ -162,6 +184,7 @@ namespace storm { operatorAtTopLevelRequired = false; multiObjectiveFormulaAtTopLevelRequired = false; operatorsAtTopLevelOfMultiObjectiveFormulasRequired = false; + quantileFormulaAtTopLevelRequired = false; rewardAccumulation = false; } @@ -215,6 +238,15 @@ namespace storm { return *this; } + bool FragmentSpecification::areQuantileFormulasAllowed() const { + return quantileFormula; + } + + FragmentSpecification& FragmentSpecification::setQuantileFormulasAllowed( bool newValue) { + this->quantileFormula = newValue; + return *this; + } + bool FragmentSpecification::areGloballyFormulasAllowed() const { return globallyFormula; } @@ -566,7 +598,16 @@ namespace storm { operatorsAtTopLevelOfMultiObjectiveFormulasRequired = newValue; return *this; } - + + bool FragmentSpecification::isQuantileFormulaAtTopLevelRequired() const { + return quantileFormulaAtTopLevelRequired; + } + + FragmentSpecification& FragmentSpecification::setQuantileFormulaAtTopLevelRequired(bool newValue) { + quantileFormulaAtTopLevelRequired = newValue; + return *this; + } + bool FragmentSpecification::isRewardAccumulationAllowed() const { return rewardAccumulation; } diff --git a/src/storm/logic/FragmentSpecification.h b/src/storm/logic/FragmentSpecification.h index 012716c78..aa48c0738 100644 --- a/src/storm/logic/FragmentSpecification.h +++ b/src/storm/logic/FragmentSpecification.h @@ -33,6 +33,9 @@ namespace storm { bool areMultiObjectiveFormulasAllowed() const; FragmentSpecification& setMultiObjectiveFormulasAllowed( bool newValue); + + bool areQuantileFormulasAllowed() const; + FragmentSpecification& setQuantileFormulasAllowed( bool newValue); bool areGloballyFormulasAllowed() const; FragmentSpecification& setGloballyFormulasAllowed(bool newValue); @@ -145,6 +148,9 @@ namespace storm { bool areOperatorsAtTopLevelOfMultiObjectiveFormulasRequired() const; FragmentSpecification& setOperatorsAtTopLevelOfMultiObjectiveFormulasRequired(bool newValue); + bool isQuantileFormulaAtTopLevelRequired() const; + FragmentSpecification& setQuantileFormulaAtTopLevelRequired(bool newValue); + bool isRewardAccumulationAllowed() const; FragmentSpecification& setRewardAccumulationAllowed(bool newValue); @@ -161,6 +167,7 @@ namespace storm { bool longRunAverageOperator; bool multiObjectiveFormula; + bool quantileFormula; bool globallyFormula; bool reachabilityProbabilityFormula; @@ -204,6 +211,7 @@ namespace storm { bool qualitativeOperatorResults; bool operatorAtTopLevelRequired; bool multiObjectiveFormulaAtTopLevelRequired; + bool quantileFormulaAtTopLevelRequired; bool operatorsAtTopLevelOfMultiObjectiveFormulasRequired; bool rewardAccumulation; @@ -232,6 +240,9 @@ namespace storm { // Multi-Objective formulas. FragmentSpecification multiObjective(); + + // Quantile formulas. + FragmentSpecification quantiles(); } } diff --git a/src/storm/logic/LiftableTransitionRewardsVisitor.cpp b/src/storm/logic/LiftableTransitionRewardsVisitor.cpp index db8df4962..99eeb444b 100644 --- a/src/storm/logic/LiftableTransitionRewardsVisitor.cpp +++ b/src/storm/logic/LiftableTransitionRewardsVisitor.cpp @@ -95,6 +95,10 @@ namespace storm { return result; } + boost::any LiftableTransitionRewardsVisitor::visit(QuantileFormula const& f, boost::any const& data) const { + return f.getSubformula().accept(*this, data); + } + boost::any LiftableTransitionRewardsVisitor::visit(NextFormula const& f, boost::any const& data) const { return boost::any_cast(f.getSubformula().accept(*this, data)); } diff --git a/src/storm/logic/LiftableTransitionRewardsVisitor.h b/src/storm/logic/LiftableTransitionRewardsVisitor.h index 9147ea3c7..c09f9c25b 100644 --- a/src/storm/logic/LiftableTransitionRewardsVisitor.h +++ b/src/storm/logic/LiftableTransitionRewardsVisitor.h @@ -32,6 +32,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const override; virtual boost::any visit(NextFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/QuantileFormula.cpp b/src/storm/logic/QuantileFormula.cpp new file mode 100644 index 000000000..ea8677afd --- /dev/null +++ b/src/storm/logic/QuantileFormula.cpp @@ -0,0 +1,96 @@ +#include "storm/logic/QuantileFormula.h" + +#include "storm/logic/FormulaVisitor.h" +#include "storm/utility/macros.h" +#include "storm/exceptions/InvalidArgumentException.h" + +namespace storm { + namespace logic { + + QuantileFormula::QuantileFormula(std::vector> const& boundVariables, std::shared_ptr subformula) : boundVariables(boundVariables), subformula(subformula) { + STORM_LOG_THROW(!boundVariables.empty(), storm::exceptions::InvalidArgumentException, "Quantile formula without bound variables are invalid."); + } + + QuantileFormula::~QuantileFormula() { + // Intentionally left empty + } + + bool QuantileFormula::isQuantileFormula() const { + return true; + } + + bool QuantileFormula::hasQuantitativeResult() const { + return true; + } + + bool QuantileFormula::hasNumericalResult() const { + return !isMultiDimensional(); + } + + bool QuantileFormula::hasParetoCurveResult() const { + return isMultiDimensional(); + } + + Formula const& QuantileFormula::getSubformula() const { + return *subformula; + } + + uint64_t QuantileFormula::getDimension() const { + return boundVariables.size(); + } + + bool QuantileFormula::isMultiDimensional() const { + return getDimension() > 1; + } + + storm::expressions::Variable const& QuantileFormula::getBoundVariable() const { + STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique bound variables. However, there are multiple bound variables defined."); + return boundVariables.front().second; + } + + storm::expressions::Variable const& QuantileFormula::getBoundVariable(uint64_t index) const { + STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested bound variable with index" << index << ". However, there are only " << boundVariables.size() << " bound variables."); + return boundVariables[index].second; + } + + std::vector> const& QuantileFormula::getBoundVariables() const { + return boundVariables; + } + + storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection() const { + STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique optimization direction of the bound variables. However, there are multiple bound variables defined."); + return boundVariables.front().first; + } + + storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection(uint64_t index) const { + STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested optimization direction with index" << index << ". However, there are only " << boundVariables.size() << " bound variables."); + return boundVariables[index].first; + } + + boost::any QuantileFormula::accept(FormulaVisitor const& visitor, boost::any const& data) const { + return visitor.visit(*this, data); + } + + void QuantileFormula::gatherAtomicExpressionFormulas(std::vector>& atomicExpressionFormulas) const { + subformula->gatherAtomicExpressionFormulas(atomicExpressionFormulas); + } + + void QuantileFormula::gatherAtomicLabelFormulas(std::vector>& atomicLabelFormulas) const { + subformula->gatherAtomicLabelFormulas(atomicLabelFormulas); + } + + void QuantileFormula::gatherReferencedRewardModels(std::set& referencedRewardModels) const { + subformula->gatherReferencedRewardModels(referencedRewardModels); + } + + std::ostream& QuantileFormula::writeToStream(std::ostream& out) const { + out << "quantile("; + for (auto const& bv : boundVariables) { + out << bv.first << " " << bv.second.getName() << ", "; + } + subformula->writeToStream(out); + out << ")"; + return out; + } + } +} diff --git a/src/storm/logic/QuantileFormula.h b/src/storm/logic/QuantileFormula.h new file mode 100644 index 000000000..2690171a5 --- /dev/null +++ b/src/storm/logic/QuantileFormula.h @@ -0,0 +1,41 @@ +#pragma once + +#include "storm/logic/Formula.h" +#include "storm/solver/OptimizationDirection.h" + +namespace storm { + namespace logic { + class QuantileFormula : public Formula { + public: + QuantileFormula(std::vector> const& boundVariables, std::shared_ptr subformula); + + virtual ~QuantileFormula(); + + virtual bool isQuantileFormula() const override; + + virtual bool hasQuantitativeResult() const override; // Result is numerical or a pareto curve + virtual bool hasNumericalResult() const; // Result is numerical + virtual bool hasParetoCurveResult() const; // Result is a pareto curve + + Formula const& getSubformula() const; + uint64_t getDimension() const; + bool isMultiDimensional() const; + + storm::expressions::Variable const& getBoundVariable() const; + storm::expressions::Variable const& getBoundVariable(uint64_t index) const; + std::vector> const& getBoundVariables() const; + storm::solver::OptimizationDirection const& getOptimizationDirection() const; + storm::solver::OptimizationDirection const& getOptimizationDirection(uint64_t index) const; + + virtual boost::any accept(FormulaVisitor const& visitor, boost::any const& data) const override; + virtual void gatherAtomicExpressionFormulas(std::vector>& atomicExpressionFormulas) const override; + virtual void gatherAtomicLabelFormulas(std::vector>& atomicLabelFormulas) const override; + virtual void gatherReferencedRewardModels(std::set& referencedRewardModels) const override; + + virtual std::ostream& writeToStream(std::ostream& out) const override; + private: + std::vector> boundVariables; + std::shared_ptr subformula; + }; + } +} diff --git a/src/storm/logic/ToExpressionVisitor.cpp b/src/storm/logic/ToExpressionVisitor.cpp index ddfe0ea6d..27ea1fa6b 100644 --- a/src/storm/logic/ToExpressionVisitor.cpp +++ b/src/storm/logic/ToExpressionVisitor.cpp @@ -87,6 +87,10 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); } + boost::any ToExpressionVisitor::visit(QuantileFormula const&, boost::any const&) const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); + } + boost::any ToExpressionVisitor::visit(NextFormula const&, boost::any const&) const { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); } diff --git a/src/storm/logic/ToExpressionVisitor.h b/src/storm/logic/ToExpressionVisitor.h index 4fe91bad3..345698f5a 100644 --- a/src/storm/logic/ToExpressionVisitor.h +++ b/src/storm/logic/ToExpressionVisitor.h @@ -26,6 +26,7 @@ namespace storm { virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageRewardFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(QuantileFormula const& f, boost::any const& data) const override; virtual boost::any visit(NextFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/storage/jani/JSONExporter.cpp b/src/storm/storage/jani/JSONExporter.cpp index 1aa2522fb..fce7562b3 100644 --- a/src/storm/storage/jani/JSONExporter.cpp +++ b/src/storm/storage/jani/JSONExporter.cpp @@ -411,6 +411,10 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Jani currently does not support conversion of a multi-objective formula"); } + boost::any FormulaToJaniJson::visit(storm::logic::QuantileFormula const&, boost::any const&) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Jani currently does not support conversion of a Quantile formula"); + } + boost::any FormulaToJaniJson::visit(storm::logic::NextFormula const& f, boost::any const& data) const { modernjson::json opDecl; opDecl["op"] = "U"; diff --git a/src/storm/storage/jani/JSONExporter.h b/src/storm/storage/jani/JSONExporter.h index 0bd1486f2..9a99a443a 100644 --- a/src/storm/storage/jani/JSONExporter.h +++ b/src/storm/storage/jani/JSONExporter.h @@ -64,6 +64,7 @@ namespace storm { virtual boost::any visit(storm::logic::LongRunAverageOperatorFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::LongRunAverageRewardFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::MultiObjectiveFormula const& f, boost::any const& data) const; + virtual boost::any visit(storm::logic::QuantileFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::NextFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::ProbabilityOperatorFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::RewardOperatorFormula const& f, boost::any const& data) const; From 9e282c8bb21e39b6fe56e494b12b2f5dd1ab664f Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 25 Jan 2019 10:53:28 +0100 Subject: [PATCH 02/20] QuantileFormulas: A boost::spiritual abyss. --- .../parser/FormulaParserGrammar.cpp | 27 ++++++++++++++++--- .../parser/FormulaParserGrammar.h | 7 ++++- 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/src/storm-parsers/parser/FormulaParserGrammar.cpp b/src/storm-parsers/parser/FormulaParserGrammar.cpp index ceca76e8e..6f565fc20 100644 --- a/src/storm-parsers/parser/FormulaParserGrammar.cpp +++ b/src/storm-parsers/parser/FormulaParserGrammar.cpp @@ -123,13 +123,18 @@ namespace storm { multiFormula = (qi::lit("multi") > qi::lit("(") >> ((pathFormula(storm::logic::FormulaContext::Probability) | stateFormula) % qi::lit(",")) >> qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createMultiFormula, phoenix::ref(*this), qi::_1)]; multiFormula.name("Multi formula"); - - stateFormula = (orStateFormula | multiFormula); - stateFormula.name("state formula"); - + identifier %= qi::as_string[qi::raw[qi::lexeme[((qi::alpha | qi::char_('_') | qi::char_('.')) >> *(qi::alnum | qi::char_('_')))]]]; identifier.name("identifier"); + quantileBoundVariable = ((qi::lit("min")[qi::_a = storm::solver::OptimizationDirection::Minimize] | qi::lit("max")[qi::_a = storm::solver::OptimizationDirection::Maximize]) >> identifier)[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileBoundVariables, phoenix::ref(*this), qi::_a, qi::_1)]; + quantileBoundVariable.name("Quantile bound variable"); + quantileFormula = (qi::lit("quantile") > qi::lit("(") >> (quantileBoundVariable % qi::lit(",")) >> qi::lit(",") >> stateFormula > qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileFormula, phoenix::ref(*this), qi::_1, qi::_2)]; + quantileFormula.name("Quantile formula"); + + stateFormula = (orStateFormula | multiFormula | quantileFormula); + stateFormula.name("state formula"); + formulaName = qi::lit("\"") >> identifier >> qi::lit("\"") >> qi::lit(":"); formulaName.name("formula name"); @@ -419,6 +424,20 @@ namespace storm { } } + std::pair FormulaParserGrammar::createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName) { + STORM_LOG_ASSERT(manager, "Mutable expression manager required to define quantile bound variable."); + STORM_LOG_THROW(!manager->hasVariable(variableName), storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in property: variable already exists."); + + storm::expressions::Variable var = manager->declareRationalVariable(variableName); + addIdentifierExpression(variableName, var); + std::pair result(dir, var); + return result; + } + + std::shared_ptr FormulaParserGrammar::createQuantileFormula(std::vector> const& boundVariables, std::shared_ptr const& subformula) { + return std::shared_ptr(new storm::logic::QuantileFormula(boundVariables, subformula)); + } + std::set FormulaParserGrammar::getUndefinedConstants(std::shared_ptr const& formula) const { std::set result; std::set usedVariables = formula->getUsedVariables(); diff --git a/src/storm-parsers/parser/FormulaParserGrammar.h b/src/storm-parsers/parser/FormulaParserGrammar.h index 75012aaf4..e582abe40 100644 --- a/src/storm-parsers/parser/FormulaParserGrammar.h +++ b/src/storm-parsers/parser/FormulaParserGrammar.h @@ -50,7 +50,8 @@ namespace storm { ("F", 5) ("G", 6) ("X", 7) - ("multi", 8); + ("multi", 8) + ("quantile", 9); } }; @@ -196,6 +197,8 @@ namespace storm { qi::rule(), Skipper> longRunAverageRewardFormula; qi::rule(), Skipper> multiFormula; + qi::rule(), qi::locals, Skipper> quantileBoundVariable; + qi::rule(), Skipper> quantileFormula; // Parser that is used to recognize doubles only (as opposed to Spirit's double_ parser). boost::spirit::qi::real_parser> strict_double; @@ -229,6 +232,8 @@ namespace storm { std::shared_ptr createBinaryBooleanStateFormula(std::shared_ptr const& leftSubformula, std::shared_ptr const& rightSubformula, storm::logic::BinaryBooleanStateFormula::OperatorType operatorType); std::shared_ptr createUnaryBooleanStateFormula(std::shared_ptr const& subformula, boost::optional const& operatorType); std::shared_ptr createMultiFormula(std::vector> const& subformulas); + std::pair createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName); + std::shared_ptr createQuantileFormula(std::vector> const& boundVariables, std::shared_ptr const& subformula); std::set getUndefinedConstants(std::shared_ptr const& formula) const; storm::jani::Property createProperty(boost::optional const& propertyName, storm::modelchecker::FilterType const& filterType, std::shared_ptr const& formula, std::shared_ptr const& states); From fe535b9ba684d89b8f0a495aa3616de58ed5ab7c Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 25 Jan 2019 10:54:56 +0100 Subject: [PATCH 03/20] Fixed QuantileFormula's writeToStream. --- src/storm/logic/QuantileFormula.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/storm/logic/QuantileFormula.cpp b/src/storm/logic/QuantileFormula.cpp index ea8677afd..322f36fdd 100644 --- a/src/storm/logic/QuantileFormula.cpp +++ b/src/storm/logic/QuantileFormula.cpp @@ -86,7 +86,7 @@ namespace storm { std::ostream& QuantileFormula::writeToStream(std::ostream& out) const { out << "quantile("; for (auto const& bv : boundVariables) { - out << bv.first << " " << bv.second.getName() << ", "; + out << (storm::solver::minimize(bv.first) ? "min" : "max") << " " << bv.second.getName() << ", "; } subformula->writeToStream(out); out << ")"; From dc2654ce609d0728fe2b9b3522048d0e1518239f Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 25 Jan 2019 21:08:24 +0100 Subject: [PATCH 04/20] Quantiles: made the SparseMdpPrctlModelChecker call the QuantileHelper for quantile formulas --- src/storm/logic/Formula.cpp | 8 ++++++ src/storm/logic/Formula.h | 3 +++ .../modelchecker/AbstractModelChecker.cpp | 7 +++++ src/storm/modelchecker/AbstractModelChecker.h | 5 +++- .../prctl/SparseMdpPrctlModelChecker.cpp | 26 ++++++++++++++++-- .../prctl/SparseMdpPrctlModelChecker.h | 1 + .../prctl/helper/rewardbounded/Dimension.h | 2 ++ .../MultiDimensionalRewardUnfolding.cpp | 27 ++++++++++--------- .../helper/rewardbounded/QuantileHelper.cpp | 10 +++++++ .../helper/rewardbounded/QuantileHelper.h | 21 +++++++++++++++ .../results/ExplicitParetoCurveCheckResult.h | 4 +-- .../results/ParetoCurveCheckResult.cpp | 22 ++++++++++++--- .../results/ParetoCurveCheckResult.h | 6 +++-- 13 files changed, 120 insertions(+), 22 deletions(-) create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h diff --git a/src/storm/logic/Formula.cpp b/src/storm/logic/Formula.cpp index fe3a1e41a..58970c77b 100644 --- a/src/storm/logic/Formula.cpp +++ b/src/storm/logic/Formula.cpp @@ -209,6 +209,14 @@ namespace storm { return dynamic_cast(*this); } + QuantileFormula& Formula::asQuantileFormula() { + return dynamic_cast(*this); + } + + QuantileFormula const& Formula::asQuantileFormula() const { + return dynamic_cast(*this); + } + BinaryStateFormula& Formula::asBinaryStateFormula() { return dynamic_cast(*this); } diff --git a/src/storm/logic/Formula.h b/src/storm/logic/Formula.h index bd88374a3..dc6396e71 100644 --- a/src/storm/logic/Formula.h +++ b/src/storm/logic/Formula.h @@ -110,6 +110,9 @@ namespace storm { MultiObjectiveFormula& asMultiObjectiveFormula(); MultiObjectiveFormula const& asMultiObjectiveFormula() const; + QuantileFormula& asQuantileFormula(); + QuantileFormula const& asQuantileFormula() const; + BinaryStateFormula& asBinaryStateFormula(); BinaryStateFormula const& asBinaryStateFormula() const; diff --git a/src/storm/modelchecker/AbstractModelChecker.cpp b/src/storm/modelchecker/AbstractModelChecker.cpp index 42bb2a753..2c12639c3 100644 --- a/src/storm/modelchecker/AbstractModelChecker.cpp +++ b/src/storm/modelchecker/AbstractModelChecker.cpp @@ -48,6 +48,8 @@ namespace storm { return this->checkStateFormula(env, checkTask.substituteFormula(formula.asStateFormula())); } else if (formula.isMultiObjectiveFormula()){ return this->checkMultiObjectiveFormula(env, checkTask.substituteFormula(formula.asMultiObjectiveFormula())); + } else if (formula.isQuantileFormula()){ + return this->checkQuantileFormula(env, checkTask.substituteFormula(formula.asQuantileFormula())); } STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << formula << "' is invalid."); } @@ -311,6 +313,11 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker (" << getClassName() << ") does not support the formula: " << checkTask.getFormula() << "."); } + template + std::unique_ptr AbstractModelChecker::checkQuantileFormula(Environment const& env, CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker (" << getClassName() << ") does not support the formula: " << checkTask.getFormula() << "."); + } + /////////////////////////////////////////////// // Explicitly instantiate the template class. /////////////////////////////////////////////// diff --git a/src/storm/modelchecker/AbstractModelChecker.h b/src/storm/modelchecker/AbstractModelChecker.h index 44373de03..80bf18fb6 100644 --- a/src/storm/modelchecker/AbstractModelChecker.h +++ b/src/storm/modelchecker/AbstractModelChecker.h @@ -91,7 +91,10 @@ namespace storm { // The methods to check multi-objective formulas. virtual std::unique_ptr checkMultiObjectiveFormula(Environment const& env, CheckTask const& checkTask); - + + // The methods to check quantile formulas. + virtual std::unique_ptr checkQuantileFormula(Environment const& env, CheckTask const& checkTask); + }; } } diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index b78ad3856..e83fba585 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -7,6 +7,7 @@ #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h" #include "storm/logic/FragmentSpecification.h" @@ -14,6 +15,7 @@ #include "storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" +#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h" #include "storm/modelchecker/multiobjective/multiObjectiveModelChecking.h" #include "storm/solver/SolveGoal.h" @@ -40,13 +42,17 @@ namespace storm { storm::logic::Formula const& formula = checkTask.getFormula(); if (formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setTotalRewardFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true).setTimeOperatorsAllowed(true).setReachbilityTimeFormulasAllowed(true))) { return true; - } else { + } else if (formula.isInFragment(storm::logic::multiObjective().setCumulativeRewardFormulasAllowed(true).setTimeBoundedCumulativeRewardFormulasAllowed(true).setStepBoundedCumulativeRewardFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setTimeBoundedUntilFormulasAllowed(true).setStepBoundedUntilFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true))) { // Check whether we consider a multi-objective formula // For multi-objective model checking, each initial state requires an individual scheduler (in contrast to single-objective model checking). Let's exclude multiple initial states. if (this->getModel().getInitialStates().getNumberOfSetBits() > 1) return false; if (!checkTask.isOnlyInitialStatesRelevantSet()) return false; - return formula.isInFragment(storm::logic::multiObjective().setCumulativeRewardFormulasAllowed(true).setTimeBoundedCumulativeRewardFormulasAllowed(true).setStepBoundedCumulativeRewardFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setTimeBoundedUntilFormulasAllowed(true).setStepBoundedUntilFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true)); + return true; + } else if (formula.isInFragment(storm::logic::quantiles())) { + if (this->getModel().getInitialStates().getNumberOfSetBits() > 1) return false; + return true; } + return false; } template @@ -220,6 +226,22 @@ namespace storm { return multiobjective::performMultiObjectiveModelChecking(env, this->getModel(), checkTask.getFormula()); } + template + std::unique_ptr SparseMdpPrctlModelChecker::checkQuantileFormula(Environment const& env, CheckTask const& checkTask) { + STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for models with a single initial states."); + STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, "Quantiles not supported on models with multiple initial states."); + uint64_t initialState = *this->getModel().getInitialStates().begin(); + + helper::rewardbounded::QuantileHelper qHelper(this->getModel(), checkTask.getFormula()); + auto res = qHelper.computeMultiDimensionalQuantile(); + + if (res.size() == 1 && res.front().size() == 1) { + return std::unique_ptr(new ExplicitQuantitativeCheckResult(initialState, std::move(res.front().front()))); + } else { + return std::unique_ptr(new ExplicitParetoCurveCheckResult(initialState, std::move(res))); + } + } + template class SparseMdpPrctlModelChecker>; #ifdef STORM_HAVE_CARL diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 19c2cc6a1..38a0a090c 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -33,6 +33,7 @@ namespace storm { virtual std::unique_ptr computeLongRunAverageProbabilities(Environment const& env, CheckTask const& checkTask) override; virtual std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr checkMultiObjectiveFormula(Environment const& env, CheckTask const& checkTask) override; + virtual std::unique_ptr checkQuantileFormula(Environment const& env, CheckTask const& checkTask) override; }; } // namespace modelchecker diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h index 80059ce1e..76127e5e8 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -3,6 +3,7 @@ #include #include "storm/storage/BitVector.h" +#include "storm/solver/OptimizationDirection.h" namespace storm { namespace modelchecker { @@ -18,6 +19,7 @@ namespace storm { ValueType scalingFactor; storm::storage::BitVector dependentDimensions; boost::optional maxValue; + boost::optional optimizationDirection; }; } } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index efc3c81f7..87d2b968d 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -88,6 +88,7 @@ namespace storm { dimension.memoryLabel = memLabel; dimension.isUpperBounded = subformula.hasUpperBound(dim); // for simplicity we do not allow intervals or unbounded formulas. + // TODO: Quantiles: allow unbounded formulas STORM_LOG_THROW(subformula.hasLowerBound(dim) != dimension.isUpperBounded, 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."); @@ -216,20 +217,22 @@ namespace storm { bound = dimFormula.asCumulativeRewardFormula().getBound(); isStrict = dimFormula.asCumulativeRewardFormula().isBoundStrict(); } - STORM_LOG_THROW(!bound.containsVariables(), storm::exceptions::NotSupportedException, "The bound " << bound << " contains undefined constants."); - ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); - STORM_LOG_THROW(dimensions[dim].isUpperBounded || isStrict || !storm::utility::isZero(discretizedBound), storm::exceptions::NotSupportedException, "Lower bounds need to be either strict or greater than zero."); - discretizedBound /= dimensions[dim].scalingFactor; - if (storm::utility::isInteger(discretizedBound)) { - if (isStrict == dimensions[dim].isUpperBounded) { - discretizedBound -= storm::utility::one(); + + if (bound.containsVariables()) { + ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); + STORM_LOG_THROW(dimensions[dim].isUpperBounded || isStrict || !storm::utility::isZero(discretizedBound), storm::exceptions::NotSupportedException, "Lower bounds need to be either strict or greater than zero."); + discretizedBound /= dimensions[dim].scalingFactor; + if (storm::utility::isInteger(discretizedBound)) { + if (isStrict == dimensions[dim].isUpperBounded) { + discretizedBound -= storm::utility::one(); + } + } else { + discretizedBound = storm::utility::floor(discretizedBound); } - } else { - discretizedBound = storm::utility::floor(discretizedBound); + uint64_t dimensionValue = storm::utility::convertNumber(discretizedBound); + STORM_LOG_THROW(epochManager.isValidDimensionValue(dimensionValue), storm::exceptions::NotSupportedException, "The bound " << bound << " is too high for the considered number of dimensions."); + dimensions[dim].maxValue = dimensionValue; } - uint64_t dimensionValue = storm::utility::convertNumber(discretizedBound); - STORM_LOG_THROW(epochManager.isValidDimensionValue(dimensionValue), storm::exceptions::NotSupportedException, "The bound " << bound << " is too high for the considered number of dimensions."); - dimensions[dim].maxValue = dimensionValue; } } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp new file mode 100644 index 000000000..a502569eb --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -0,0 +1,10 @@ +#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h" + +namespace storm { + namespace modelchecker { + namespace helper { + namespace rewardbounded { + } + } + } +} diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h new file mode 100644 index 000000000..4d31b1624 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -0,0 +1,21 @@ +#pragma once + +#include "storm/logic/QuantileFormula.h" + +namespace storm { + namespace modelchecker { + namespace helper { + namespace rewardbounded { + + template + class QuantileHelper { + typedef typename ModelType::ValueType ValueType; + public: + QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& formula) {} + + std::vector> computeMultiDimensionalQuantile() { return {{27}};} + }; + } + } + } +} diff --git a/src/storm/modelchecker/results/ExplicitParetoCurveCheckResult.h b/src/storm/modelchecker/results/ExplicitParetoCurveCheckResult.h index 5bc42feff..e25221958 100644 --- a/src/storm/modelchecker/results/ExplicitParetoCurveCheckResult.h +++ b/src/storm/modelchecker/results/ExplicitParetoCurveCheckResult.h @@ -12,8 +12,8 @@ namespace storm { class ExplicitParetoCurveCheckResult : public ParetoCurveCheckResult { public: ExplicitParetoCurveCheckResult(); - ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector::point_type> const& points, typename ParetoCurveCheckResult::polytope_type const& underApproximation, typename ParetoCurveCheckResult::polytope_type const& overApproximation); - ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector::point_type>&& points, typename ParetoCurveCheckResult::polytope_type&& underApproximation, typename ParetoCurveCheckResult::polytope_type&& overApproximation); + ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector::point_type> const& points, typename ParetoCurveCheckResult::polytope_type const& underApproximation = nullptr, typename ParetoCurveCheckResult::polytope_type const& overApproximation = nullptr); + ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector::point_type>&& points, typename ParetoCurveCheckResult::polytope_type&& underApproximation = nullptr, typename ParetoCurveCheckResult::polytope_type&& overApproximation = nullptr); ExplicitParetoCurveCheckResult(ExplicitParetoCurveCheckResult const& other) = default; ExplicitParetoCurveCheckResult& operator=(ExplicitParetoCurveCheckResult const& other) = default; diff --git a/src/storm/modelchecker/results/ParetoCurveCheckResult.cpp b/src/storm/modelchecker/results/ParetoCurveCheckResult.cpp index 39ae4e265..ffc34528a 100644 --- a/src/storm/modelchecker/results/ParetoCurveCheckResult.cpp +++ b/src/storm/modelchecker/results/ParetoCurveCheckResult.cpp @@ -31,22 +31,38 @@ namespace storm { return points; } + template + bool ParetoCurveCheckResult::hasUnderApproximation() const { + return bool(underApproximation); + } + + template + bool ParetoCurveCheckResult::hasOverApproximation() const { + return bool(overApproximation); + } + template typename ParetoCurveCheckResult::polytope_type const& ParetoCurveCheckResult::getUnderApproximation() const { + STORM_LOG_ASSERT(hasUnderApproximation(), "Requested under approx. of Pareto curve although it does not exist."); return underApproximation; } template typename ParetoCurveCheckResult::polytope_type const& ParetoCurveCheckResult::getOverApproximation() const { + STORM_LOG_ASSERT(hasUnderApproximation(), "Requested over approx. of Pareto curve although it does not exist."); return overApproximation; } template std::ostream& ParetoCurveCheckResult::writeToStream(std::ostream& out) const { out << std::endl; - out << "Underapproximation of achievable values: " << underApproximation->toString() << std::endl; - out << "Overapproximation of achievable values: " << overApproximation->toString() << std::endl; - out << points.size() << " pareto optimal points found (Note that these points are safe, i.e., contained in the underapproximation, but there is no guarantee for optimality):" << std::endl; + if (hasUnderApproximation()) { + out << "Underapproximation of achievable values: " << underApproximation->toString() << std::endl; + } + if (hasOverApproximation()) { + out << "Overapproximation of achievable values: " << overApproximation->toString() << std::endl; + } + out << points.size() << " Pareto optimal points found:" << std::endl; for(auto const& p : points) { out << " ("; for(auto it = p.begin(); it != p.end(); ++it){ diff --git a/src/storm/modelchecker/results/ParetoCurveCheckResult.h b/src/storm/modelchecker/results/ParetoCurveCheckResult.h index e2b3be825..9397bc4c5 100644 --- a/src/storm/modelchecker/results/ParetoCurveCheckResult.h +++ b/src/storm/modelchecker/results/ParetoCurveCheckResult.h @@ -19,14 +19,16 @@ namespace storm { virtual bool isParetoCurveCheckResult() const override; std::vector const& getPoints() const; + bool hasUnderApproximation() const; + bool hasOverApproximation() const; polytope_type const& getUnderApproximation() const; polytope_type const& getOverApproximation() const; virtual std::ostream& writeToStream(std::ostream& out) const override; protected: - ParetoCurveCheckResult(std::vector const& points, polytope_type const& underApproximation, polytope_type const& overApproximation); - ParetoCurveCheckResult(std::vector&& points, polytope_type&& underApproximation, polytope_type&& overApproximation); + ParetoCurveCheckResult(std::vector const& points, polytope_type const& underApproximation = nullptr, polytope_type const& overApproximation = nullptr); + ParetoCurveCheckResult(std::vector&& points, polytope_type&& underApproximation = nullptr, polytope_type&& overApproximation = nullptr); // The pareto optimal points that have been found. std::vector points; From d3abeb5f45d74f3d9b0e6fef79fcf420e4824903 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Mon, 28 Jan 2019 12:58:18 +0100 Subject: [PATCH 05/20] Started implementation on quantiles. --- src/storm/logic/TimeBound.h | 3 +- src/storm/logic/TimeBoundType.h | 2 + .../prctl/SparseMdpPrctlModelChecker.cpp | 2 +- .../helper/rewardbounded/QuantileHelper.cpp | 81 +++++++++++++++++++ .../helper/rewardbounded/QuantileHelper.h | 16 +++- 5 files changed, 99 insertions(+), 5 deletions(-) diff --git a/src/storm/logic/TimeBound.h b/src/storm/logic/TimeBound.h index 17a4c2b41..6ab0c29c8 100644 --- a/src/storm/logic/TimeBound.h +++ b/src/storm/logic/TimeBound.h @@ -8,7 +8,8 @@ namespace storm { class TimeBound { public: TimeBound(bool strict, storm::expressions::Expression const& bound); - + TimeBound(TimeBound const& other) = default; + storm::expressions::Expression const& getBound() const; bool isStrict() const; diff --git a/src/storm/logic/TimeBoundType.h b/src/storm/logic/TimeBoundType.h index 149f76215..91f4b134c 100644 --- a/src/storm/logic/TimeBoundType.h +++ b/src/storm/logic/TimeBoundType.h @@ -29,6 +29,8 @@ namespace storm { assert(rewardName != ""); // Empty reward name is reserved. } + TimeBoundReference(TimeBoundReference const& other) = default; + TimeBoundType getType() const { return type; } diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index e83fba585..b16bd7b3c 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -233,7 +233,7 @@ namespace storm { uint64_t initialState = *this->getModel().getInitialStates().begin(); helper::rewardbounded::QuantileHelper qHelper(this->getModel(), checkTask.getFormula()); - auto res = qHelper.computeMultiDimensionalQuantile(); + auto res = qHelper.computeMultiDimensionalQuantile(env); if (res.size() == 1 && res.front().size() == 1) { return std::unique_ptr(new ExplicitQuantitativeCheckResult(initialState, std::move(res.front().front()))); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index a502569eb..54ca1622f 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -1,9 +1,90 @@ #include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h" +#include +#include +#include +#include + +#include "storm/models/sparse/Mdp.h" +#include "storm/storage/expressions/ExpressionManager.h" + +#include "storm/logic/ProbabilityOperatorFormula.h" +#include "storm/logic/BoundedUntilFormula.h" + namespace storm { namespace modelchecker { namespace helper { namespace rewardbounded { + + template + QuantileHelper::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) { + // Intentionally left empty + } + + std::shared_ptr replaceBoundsByGreaterEqZero(std::shared_ptr const& boundedUntilOperator, std::set dimensionsToReplace) { + auto const& origBoundedUntil = boundedUntilOperator->getSubformula().asBoundedUntilFormula(); + STORM_LOG_ASSERT(*(--dimensionsToReplace.end()) < 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; + 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()); + timeBoundReferences.push_back(origBoundedUntil.getTimeBoundReference(dim)); + if (dimensionsToReplace.count(dim) == 0) { + if (origBoundedUntil.hasLowerBound()) { + lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); + } else { + lowerBounds.push_back(boost::none); + } + if (origBoundedUntil.hasUpperBound()) { + upperBounds.emplace_back(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim)); + } else { + upperBounds.push_back(boost::none); + } + } else { + storm::expressions::Expression zero; + if (origBoundedUntil.hasLowerBound(dim)) { + zero = origBoundedUntil.getLowerBound(dim).getManager().rational(0.0); + } else { + STORM_LOG_THROW(origBoundedUntil.hasUpperBound(dim), storm::exceptions::InvalidOperationException, "The given bounded until formula has no cost-bound for one dimension."); + zero = origBoundedUntil.getUpperBound(dim).getManager().rational(0.0); + } + lowerBounds.emplace_back(false, zero); + upperBounds.push_back(boost::none); + } + } + auto newBoundedUntil = std::make_shared(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences); + return std::make_shared(newBoundedUntil, boundedUntilOperator->getOperatorInformation()); + } + + + + template + std::vector> QuantileHelper::computeMultiDimensionalQuantile(Environment const& env) { + std::vector> result; + if (!computeUnboundedValue(env)) { + STORM_LOG_INFO("Bound is not satisfiable."); + result.emplace_back(); + for (auto const& bv : quantileFormula.getBoundVariables()) { + result.front().push_back(storm::solver::minimize(bv.first) ? storm::utility::infinity() : -storm::utility::infinity()); + } + } else { + result = {{27}}; + } + return result; + } + + template + bool QuantileHelper::computeUnboundedValue(Environment const& env) { + return false; + } + + + 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 4d31b1624..276b9faed 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -3,6 +3,8 @@ #include "storm/logic/QuantileFormula.h" namespace storm { + class Environment; + namespace modelchecker { namespace helper { namespace rewardbounded { @@ -11,9 +13,17 @@ namespace storm { class QuantileHelper { typedef typename ModelType::ValueType ValueType; public: - QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& formula) {} - - std::vector> computeMultiDimensionalQuantile() { return {{27}};} + QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula); + + std::vector> computeMultiDimensionalQuantile(Environment const& env); + + + + private: + bool computeUnboundedValue(Environment const& env); + + ModelType const& model; + storm::logic::QuantileFormula const& quantileFormula; }; } } From d796ee74de868e31bfdcd96775e42599fca368b4 Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 29 Jan 2019 10:28:59 +0100 Subject: [PATCH 06/20] (workplace switch) --- .../prctl/helper/rewardbounded/QuantileHelper.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 54ca1622f..336c1d299 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -21,8 +21,8 @@ namespace storm { // Intentionally left empty } - std::shared_ptr replaceBoundsByGreaterEqZero(std::shared_ptr const& boundedUntilOperator, std::set dimensionsToReplace) { - auto const& origBoundedUntil = boundedUntilOperator->getSubformula().asBoundedUntilFormula(); + std::shared_ptr replaceBoundsByGreaterEqZero(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::set const& dimensionsToReplace) { + auto const& origBoundedUntil = boundedUntilOperator.getSubformula().asBoundedUntilFormula(); STORM_LOG_ASSERT(*(--dimensionsToReplace.end()) < 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; std::vector> lowerBounds, upperBounds; @@ -39,7 +39,7 @@ namespace storm { lowerBounds.push_back(boost::none); } if (origBoundedUntil.hasUpperBound()) { - upperBounds.emplace_back(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim)); + lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim))); } else { upperBounds.push_back(boost::none); } @@ -51,12 +51,12 @@ namespace storm { STORM_LOG_THROW(origBoundedUntil.hasUpperBound(dim), storm::exceptions::InvalidOperationException, "The given bounded until formula has no cost-bound for one dimension."); zero = origBoundedUntil.getUpperBound(dim).getManager().rational(0.0); } - lowerBounds.emplace_back(false, zero); + lowerBounds.push_back(storm::logic::TimeBound(false, zero)); upperBounds.push_back(boost::none); } } auto newBoundedUntil = std::make_shared(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences); - return std::make_shared(newBoundedUntil, boundedUntilOperator->getOperatorInformation()); + return std::make_shared(newBoundedUntil, boundedUntilOperator.getOperatorInformation()); } @@ -78,6 +78,7 @@ namespace storm { template bool QuantileHelper::computeUnboundedValue(Environment const& env) { + auto unboundedFormula = replaceBoundsByGreaterEqZero(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::set()); return false; } From 661d922d2e68aaf720173862e4ee79741a252fee Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 29 Jan 2019 12:56:29 +0100 Subject: [PATCH 07/20] logic/bound: Added method to compare a bound with a value. --- src/storm/logic/Bound.h | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/storm/logic/Bound.h b/src/storm/logic/Bound.h index 21f2c7282..a39edb5dc 100644 --- a/src/storm/logic/Bound.h +++ b/src/storm/logic/Bound.h @@ -16,6 +16,21 @@ namespace storm { ComparisonType comparisonType; storm::expressions::Expression threshold; + template + bool isSatisfied(ValueType const& compareValue) { + ValueType thresholdAsValueType = storm::utility::convertNumber(threshold.evaluateAsRational()); + switch(comparisonType) { + case ComparisonType::Greater: + return compareValue > thresholdAsValueType; + case ComparisonType::GreaterEqual: + return compareValue >= thresholdAsValueType; + case ComparisonType::Less: + return compareValue < thresholdAsValueType; + case ComparisonType::LessEqual: + return compareValue <= thresholdAsValueType; + } + } + friend std::ostream& operator<<(std::ostream& out, Bound const& bound); }; From 82402ba3ae301b9e16b32fde661835123aa69526 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 29 Jan 2019 13:08:20 +0100 Subject: [PATCH 08/20] rewardbounded: Moved epoch model analysis to a separate file. --- ...ewardBoundedMdpPcaaWeightVectorChecker.cpp | 2 +- .../RewardBoundedMdpPcaaWeightVectorChecker.h | 2 +- .../prctl/helper/SparseDtmcPrctlHelper.cpp | 79 +----- .../prctl/helper/SparseMdpPrctlHelper.cpp | 105 +------- .../prctl/helper/rewardbounded/EpochModel.cpp | 248 ++++++++++++++++++ .../prctl/helper/rewardbounded/EpochModel.h | 48 ++++ .../MultiDimensionalRewardUnfolding.cpp | 8 +- .../MultiDimensionalRewardUnfolding.h | 24 +- 8 files changed, 312 insertions(+), 204 deletions(-) create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h diff --git a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp index f5137b0a7..0df580eab 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp +++ b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp @@ -309,7 +309,7 @@ namespace storm { } template - void RewardBoundedMdpPcaaWeightVectorChecker::updateCachedData(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector) { + void RewardBoundedMdpPcaaWeightVectorChecker::updateCachedData(Environment const& env, helper::rewardbounded::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector) { if (epochModel.epochMatrixChanged) { // Update the cached MinMaxSolver data diff --git a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.h b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.h index 0013cdc3d..68fd403f0 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.h +++ b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.h @@ -64,7 +64,7 @@ namespace storm { void computeEpochSolution(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding::Epoch const& epoch, std::vector const& weightVector, EpochCheckingData& cachedData); - void updateCachedData(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector); + void updateCachedData(Environment const& env, typename helper::rewardbounded::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector); storm::utility::Stopwatch swAll, swEpochModelBuild, swEpochModelAnalysis; uint64_t numCheckedEpochs, numChecks; diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 0e5cf1849..28e799691 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -81,75 +81,6 @@ namespace storm { return result; } - template - std::vector analyzeTrivialDtmcEpochModel(typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { - - std::vector epochResult; - epochResult.reserve(epochModel.epochInStates.getNumberOfSetBits()); - auto stepSolutionIt = epochModel.stepSolutions.begin(); - auto stepChoiceIt = epochModel.stepChoices.begin(); - for (auto const& state : epochModel.epochInStates) { - while (*stepChoiceIt < state) { - ++stepChoiceIt; - ++stepSolutionIt; - } - if (epochModel.objectiveRewardFilter.front().get(state)) { - if (*stepChoiceIt == state) { - epochResult.push_back(epochModel.objectiveRewards.front()[state] + *stepSolutionIt); - } else { - epochResult.push_back(epochModel.objectiveRewards.front()[state]); - } - } else { - if (*stepChoiceIt == state) { - epochResult.push_back(*stepSolutionIt); - } else { - epochResult.push_back(storm::utility::zero()); - } - } - } - return epochResult; - } - - template - std::vector analyzeNonTrivialDtmcEpochModel(Environment const& env, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& linEqSolver, boost::optional const& lowerBound, boost::optional const& upperBound) { - - // Update some data for the case that the Matrix has changed - if (epochModel.epochMatrixChanged) { - x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero()); - storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; - linEqSolver = linearEquationSolverFactory.create(env, epochModel.epochMatrix); - linEqSolver->setCachingEnabled(true); - auto req = linEqSolver->getRequirements(env); - if (lowerBound) { - linEqSolver->setLowerBound(lowerBound.get()); - req.clearLowerBounds(); - } - if (upperBound) { - linEqSolver->setUpperBound(upperBound.get()); - req.clearUpperBounds(); - } - STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); - } - - // Prepare the right hand side of the equation system - b.assign(epochModel.epochMatrix.getRowCount(), storm::utility::zero()); - std::vector const& objectiveValues = epochModel.objectiveRewards.front(); - for (auto const& choice : epochModel.objectiveRewardFilter.front()) { - b[choice] = objectiveValues[choice]; - } - auto stepSolutionIt = epochModel.stepSolutions.begin(); - for (auto const& choice : epochModel.stepChoices) { - b[choice] += *stepSolutionIt; - ++stepSolutionIt; - } - assert(stepSolutionIt == epochModel.stepSolutions.end()); - - // Solve the minMax equation system - linEqSolver->solveEquations(env, x, b); - - return storm::utility::vector::filterVector(x, epochModel.epochInStates); - } - template<> std::map SparseDtmcPrctlHelper::computeRewardBoundedValues(Environment const& env, storm::models::sparse::Dtmc const& model, std::shared_ptr rewardBoundedFormula) { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The specified property is not supported by this value type."); @@ -184,8 +115,7 @@ namespace storm { // Set the correct equation problem format. storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; rewardUnfolding.setEquationSystemFormatForEpochModel(linearEquationSolverFactory.getEquationProblemFormat(preciseEnv)); - bool convertToEquationSystem = linearEquationSolverFactory.getEquationProblemFormat(preciseEnv) == solver::LinearEquationSolverProblemFormat::EquationSystem; - + storm::utility::ProgressMeasurement progress("epochs"); progress.setMaxCount(epochOrder.size()); progress.startNewMeasurement(0); @@ -194,12 +124,7 @@ namespace storm { swBuild.start(); auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); swBuild.stop(); swCheck.start(); - // If the epoch matrix is empty we do not need to solve a linear equation system - if ((convertToEquationSystem && epochModel.epochMatrix.isIdentityMatrix()) || (!convertToEquationSystem && epochModel.epochMatrix.getEntryCount() == 0)) { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialDtmcEpochModel(epochModel)); - } else { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialDtmcEpochModel(preciseEnv, epochModel, x, b, linEqSolver, lowerBound, upperBound)); - } + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(preciseEnv, x, b, linEqSolver, lowerBound, upperBound)); swCheck.stop(); if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { std::vector cdfEntry; diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 0a0ea992c..b28bcd6f0 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -90,103 +90,7 @@ namespace storm { return result; } - template - std::vector analyzeTrivialMdpEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { - // Assert that the epoch model is indeed trivial - assert(epochModel.epochMatrix.getEntryCount() == 0); - - std::vector epochResult; - epochResult.reserve(epochModel.epochInStates.getNumberOfSetBits()); - - auto stepSolutionIt = epochModel.stepSolutions.begin(); - auto stepChoiceIt = epochModel.stepChoices.begin(); - for (auto const& state : epochModel.epochInStates) { - // Obtain the best choice for this state - ValueType bestValue; - uint64_t lastChoice = epochModel.epochMatrix.getRowGroupIndices()[state + 1]; - bool isFirstChoice = true; - for (uint64_t choice = epochModel.epochMatrix.getRowGroupIndices()[state]; choice < lastChoice; ++choice) { - while (*stepChoiceIt < choice) { - ++stepChoiceIt; - ++stepSolutionIt; - } - - ValueType choiceValue = storm::utility::zero(); - if (epochModel.objectiveRewardFilter.front().get(choice)) { - choiceValue += epochModel.objectiveRewards.front()[choice]; - } - if (*stepChoiceIt == choice) { - choiceValue += *stepSolutionIt; - } - - if (isFirstChoice) { - bestValue = std::move(choiceValue); - isFirstChoice = false; - } else { - if (storm::solver::minimize(dir)) { - if (choiceValue < bestValue) { - bestValue = std::move(choiceValue); - } - } else { - if (choiceValue > bestValue) { - bestValue = std::move(choiceValue); - } - } - } - } - // Insert the solution w.r.t. this choice - epochResult.push_back(std::move(bestValue)); - } - return epochResult; - } - - template - std::vector analyzeNonTrivialMdpEpochModel(Environment const& env, OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, boost::optional const& lowerBound, boost::optional const& upperBound) { - - // Update some data for the case that the Matrix has changed - if (epochModel.epochMatrixChanged) { - x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero()); - storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxLinearEquationSolverFactory; - minMaxSolver = minMaxLinearEquationSolverFactory.create(env, epochModel.epochMatrix); - minMaxSolver->setHasUniqueSolution(); - minMaxSolver->setOptimizationDirection(dir); - minMaxSolver->setCachingEnabled(true); - minMaxSolver->setTrackScheduler(true); - auto req = minMaxSolver->getRequirements(env, dir, false); - if (lowerBound) { - minMaxSolver->setLowerBound(lowerBound.get()); - req.clearLowerBounds(); - } - if (upperBound) { - minMaxSolver->setUpperBound(upperBound.get()); - req.clearUpperBounds(); - } - STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); - minMaxSolver->setRequirementsChecked(); - } else { - auto choicesTmp = minMaxSolver->getSchedulerChoices(); - minMaxSolver->setInitialScheduler(std::move(choicesTmp)); - } - - // Prepare the right hand side of the equation system - b.assign(epochModel.epochMatrix.getRowCount(), storm::utility::zero()); - std::vector const& objectiveValues = epochModel.objectiveRewards.front(); - for (auto const& choice : epochModel.objectiveRewardFilter.front()) { - b[choice] = objectiveValues[choice]; - } - auto stepSolutionIt = epochModel.stepSolutions.begin(); - for (auto const& choice : epochModel.stepChoices) { - b[choice] += *stepSolutionIt; - ++stepSolutionIt; - } - assert(stepSolutionIt == epochModel.stepSolutions.end()); - - // Solve the minMax equation system - minMaxSolver->solveEquations(env, x, b); - - return storm::utility::vector::filterVector(x, epochModel.epochInStates); - } - + template std::map SparseMdpPrctlHelper::computeRewardBoundedValues(Environment const& env, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates) { storm::utility::Stopwatch swAll(true), swBuild, swCheck; @@ -218,12 +122,7 @@ namespace storm { swBuild.start(); auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); swBuild.stop(); swCheck.start(); - // If the epoch matrix is empty we do not need to solve a linear equation system - if (epochModel.epochMatrix.getEntryCount() == 0) { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialMdpEpochModel(dir, epochModel)); - } else { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialMdpEpochModel(preciseEnv, dir, epochModel, x, b, minMaxSolver, lowerBound, upperBound)); - } + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(preciseEnv, dir, x, b, minMaxSolver, lowerBound, upperBound)); swCheck.stop(); if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { std::vector cdfEntry; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp new file mode 100644 index 000000000..37c39e691 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp @@ -0,0 +1,248 @@ +#include "storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h" +#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" + +#include "storm/exceptions/UncheckedRequirementException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + namespace rewardbounded { + + + + template + std::vector analyzeTrivialDtmcEpochModel(EpochModel& epochModel) { + + std::vector epochResult; + epochResult.reserve(epochModel.epochInStates.getNumberOfSetBits()); + auto stepSolutionIt = epochModel.stepSolutions.begin(); + auto stepChoiceIt = epochModel.stepChoices.begin(); + for (auto const& state : epochModel.epochInStates) { + while (*stepChoiceIt < state) { + ++stepChoiceIt; + ++stepSolutionIt; + } + if (epochModel.objectiveRewardFilter.front().get(state)) { + if (*stepChoiceIt == state) { + epochResult.push_back(epochModel.objectiveRewards.front()[state] + *stepSolutionIt); + } else { + epochResult.push_back(epochModel.objectiveRewards.front()[state]); + } + } else { + if (*stepChoiceIt == state) { + epochResult.push_back(*stepSolutionIt); + } else { + epochResult.push_back(storm::utility::zero()); + } + } + } + return epochResult; + } + + template + std::vector analyzeNonTrivialDtmcEpochModel(Environment const& env, EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& linEqSolver, boost::optional const& lowerBound, boost::optional const& upperBound) { + + // Update some data for the case that the Matrix has changed + if (epochModel.epochMatrixChanged) { + x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero()); + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + linEqSolver = linearEquationSolverFactory.create(env, epochModel.epochMatrix); + linEqSolver->setCachingEnabled(true); + auto req = linEqSolver->getRequirements(env); + if (lowerBound) { + linEqSolver->setLowerBound(lowerBound.get()); + req.clearLowerBounds(); + } + if (upperBound) { + linEqSolver->setUpperBound(upperBound.get()); + req.clearUpperBounds(); + } + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); + } + + // Prepare the right hand side of the equation system + b.assign(epochModel.epochMatrix.getRowCount(), storm::utility::zero()); + std::vector const& objectiveValues = epochModel.objectiveRewards.front(); + for (auto const& choice : epochModel.objectiveRewardFilter.front()) { + b[choice] = objectiveValues[choice]; + } + auto stepSolutionIt = epochModel.stepSolutions.begin(); + for (auto const& choice : epochModel.stepChoices) { + b[choice] += *stepSolutionIt; + ++stepSolutionIt; + } + assert(stepSolutionIt == epochModel.stepSolutions.end()); + + // Solve the minMax equation system + linEqSolver->solveEquations(env, x, b); + + return storm::utility::vector::filterVector(x, epochModel.epochInStates); + } + + + + template + std::vector analyzeTrivialMdpEpochModel(OptimizationDirection dir, EpochModel& epochModel) { + // Assert that the epoch model is indeed trivial + assert(epochModel.epochMatrix.getEntryCount() == 0); + + std::vector epochResult; + epochResult.reserve(epochModel.epochInStates.getNumberOfSetBits()); + + auto stepSolutionIt = epochModel.stepSolutions.begin(); + auto stepChoiceIt = epochModel.stepChoices.begin(); + for (auto const& state : epochModel.epochInStates) { + // Obtain the best choice for this state + ValueType bestValue; + uint64_t lastChoice = epochModel.epochMatrix.getRowGroupIndices()[state + 1]; + bool isFirstChoice = true; + for (uint64_t choice = epochModel.epochMatrix.getRowGroupIndices()[state]; choice < lastChoice; ++choice) { + while (*stepChoiceIt < choice) { + ++stepChoiceIt; + ++stepSolutionIt; + } + + ValueType choiceValue = storm::utility::zero(); + if (epochModel.objectiveRewardFilter.front().get(choice)) { + choiceValue += epochModel.objectiveRewards.front()[choice]; + } + if (*stepChoiceIt == choice) { + choiceValue += *stepSolutionIt; + } + + if (isFirstChoice) { + bestValue = std::move(choiceValue); + isFirstChoice = false; + } else { + if (storm::solver::minimize(dir)) { + if (choiceValue < bestValue) { + bestValue = std::move(choiceValue); + } + } else { + if (choiceValue > bestValue) { + bestValue = std::move(choiceValue); + } + } + } + } + // Insert the solution w.r.t. this choice + epochResult.push_back(std::move(bestValue)); + } + return epochResult; + } + + template + std::vector analyzeNonTrivialMdpEpochModel(Environment const& env, OptimizationDirection dir, EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, boost::optional const& lowerBound, boost::optional const& upperBound) { + + // Update some data for the case that the Matrix has changed + if (epochModel.epochMatrixChanged) { + x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero()); + storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxLinearEquationSolverFactory; + minMaxSolver = minMaxLinearEquationSolverFactory.create(env, epochModel.epochMatrix); + minMaxSolver->setHasUniqueSolution(); + minMaxSolver->setOptimizationDirection(dir); + minMaxSolver->setCachingEnabled(true); + minMaxSolver->setTrackScheduler(true); + auto req = minMaxSolver->getRequirements(env, dir, false); + if (lowerBound) { + minMaxSolver->setLowerBound(lowerBound.get()); + req.clearLowerBounds(); + } + if (upperBound) { + minMaxSolver->setUpperBound(upperBound.get()); + req.clearUpperBounds(); + } + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); + minMaxSolver->setRequirementsChecked(); + } else { + auto choicesTmp = minMaxSolver->getSchedulerChoices(); + minMaxSolver->setInitialScheduler(std::move(choicesTmp)); + } + + // Prepare the right hand side of the equation system + b.assign(epochModel.epochMatrix.getRowCount(), storm::utility::zero()); + std::vector const& objectiveValues = epochModel.objectiveRewards.front(); + for (auto const& choice : epochModel.objectiveRewardFilter.front()) { + b[choice] = objectiveValues[choice]; + } + auto stepSolutionIt = epochModel.stepSolutions.begin(); + for (auto const& choice : epochModel.stepChoices) { + b[choice] += *stepSolutionIt; + ++stepSolutionIt; + } + assert(stepSolutionIt == epochModel.stepSolutions.end()); + + // Solve the minMax equation system + minMaxSolver->solveEquations(env, x, b); + + return storm::utility::vector::filterVector(x, epochModel.epochInStates); + } + + template<> + std::vector EpochModel::analyzeSingleObjective( + const storm::Environment &env, std::vector &x, std::vector &b, + std::unique_ptr> &linEqSolver, + const boost::optional &lowerBound, const boost::optional &upperBound) { + STORM_LOG_ASSERT(epochMatrix.hasTrivialRowGrouping(), "This operation is only allowed if no nondeterminism is present."); + STORM_LOG_ASSERT(equationSolverProblemFormat.is_initialized(), "Unknown equation problem format."); + // If the epoch matrix is empty we do not need to solve a linear equation system + bool convertToEquationSystem = (equationSolverProblemFormat == storm::solver::LinearEquationSolverProblemFormat::EquationSystem); + if ((convertToEquationSystem && epochMatrix.isIdentityMatrix()) || (!convertToEquationSystem && epochMatrix.getEntryCount() == 0)) { + return analyzeTrivialDtmcEpochModel(*this); + } else { + return analyzeNonTrivialDtmcEpochModel(env, *this, x, b, linEqSolver, lowerBound, upperBound); + } + } + + template<> + std::vector EpochModel::analyzeSingleObjective( + const storm::Environment &env, storm::OptimizationDirection dir, std::vector &x, + std::vector &b, + std::unique_ptr> &minMaxSolver, + const boost::optional &lowerBound, const boost::optional &upperBound) { + // If the epoch matrix is empty we do not need to solve a linear equation system + if (epochMatrix.getEntryCount() == 0) { + return analyzeTrivialMdpEpochModel(dir, *this); + } else { + return analyzeNonTrivialMdpEpochModel(env, dir, *this, x, b, minMaxSolver, lowerBound, upperBound); + } + } + + template<> + std::vector EpochModel::analyzeSingleObjective( + const storm::Environment &env, std::vector &x, std::vector &b, + std::unique_ptr> &linEqSolver, + const boost::optional &lowerBound, const boost::optional &upperBound) { + STORM_LOG_ASSERT(epochMatrix.hasTrivialRowGrouping(), "This operation is only allowed if no nondeterminism is present."); + STORM_LOG_ASSERT(equationSolverProblemFormat.is_initialized(), "Unknown equation problem format."); + // If the epoch matrix is empty we do not need to solve a linear equation system + bool convertToEquationSystem = (equationSolverProblemFormat == storm::solver::LinearEquationSolverProblemFormat::EquationSystem); + if ((convertToEquationSystem && epochMatrix.isIdentityMatrix()) || (!convertToEquationSystem && epochMatrix.getEntryCount() == 0)) { + return analyzeTrivialDtmcEpochModel(*this); + } else { + return analyzeNonTrivialDtmcEpochModel(env, *this, x, b, linEqSolver, lowerBound, upperBound); + } + } + + template<> + std::vector EpochModel::analyzeSingleObjective( + const storm::Environment &env, storm::OptimizationDirection dir, std::vector &x, + std::vector &b, + std::unique_ptr> &minMaxSolver, + const boost::optional &lowerBound, const boost::optional &upperBound) { + // If the epoch matrix is empty we do not need to solve a linear equation system + if (epochMatrix.getEntryCount() == 0) { + return analyzeTrivialMdpEpochModel(dir, *this); + } else { + return analyzeNonTrivialMdpEpochModel(env, dir, *this, x, b, minMaxSolver, lowerBound, upperBound); + } + } + + template class EpochModel; + template class EpochModel; + template class EpochModel; + template class EpochModel; + } + } + } +} diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h new file mode 100644 index 000000000..dae640ff5 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h @@ -0,0 +1,48 @@ +#pragma once + +#include +#include "storm/storage/SparseMatrix.h" +#include "storm/storage/BitVector.h" +#include "storm/solver/LinearEquationSolverProblemFormat.h" +#include "storm/solver/OptimizationDirection.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/solver/LinearEquationSolver.h" + +namespace storm { + class Environment; + + namespace modelchecker { + namespace helper { + namespace rewardbounded { + template + struct EpochModel { + typedef typename std::conditional>::type SolutionType; + + bool epochMatrixChanged; + storm::storage::SparseMatrix epochMatrix; + storm::storage::BitVector stepChoices; + std::vector stepSolutions; + std::vector> objectiveRewards; + std::vector objectiveRewardFilter; + storm::storage::BitVector epochInStates; + /// In case of DTMCs we have different options for the equation problem format the epoch model will have. + boost::optional equationSolverProblemFormat; + + /*! + * Analyzes the epoch model, i.e., solves the represented equation system. This method assumes a nondeterministic model. + */ + std::vector analyzeSingleObjective(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, boost::optional const& lowerBound, boost::optional const& upperBound); + + /*! + * Analyzes the epoch model, i.e., solves the represented equation system. This method assumes a deterministic model. + */ + std::vector analyzeSingleObjective(Environment const& env, std::vector& x, std::vector& b, std::unique_ptr>& linEqSolver, boost::optional const& lowerBound, boost::optional const& upperBound); + }; + + + + + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 87d2b968d..4bd08a021 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -287,7 +287,7 @@ namespace storm { } template - typename MultiDimensionalRewardUnfolding::EpochModel& MultiDimensionalRewardUnfolding::setCurrentEpoch(Epoch const& epoch) { + EpochModel& MultiDimensionalRewardUnfolding::setCurrentEpoch(Epoch const& epoch) { STORM_LOG_DEBUG("Setting model for epoch " << epochManager.toString(epoch)); // Check if we need to update the current epoch class @@ -446,8 +446,8 @@ namespace storm { if (model.isOfType(storm::models::ModelType::Dtmc)) { assert(zeroObjRewardChoices.size() == productModel->getProduct().getNumberOfStates()); assert(stepChoices.size() == productModel->getProduct().getNumberOfStates()); - STORM_LOG_ASSERT(equationSolverProblemFormatForEpochModel.is_initialized(), "Linear equation problem format was not set."); - bool convertToEquationSystem = equationSolverProblemFormatForEpochModel.get() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + STORM_LOG_ASSERT(epochModel.equationSolverProblemFormat.is_initialized(), "Linear equation problem format was not set."); + bool convertToEquationSystem = epochModel.equationSolverProblemFormat.get() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; // For DTMCs we consider the subsystem induced by the considered states. // The transitions for states with zero reward are filtered out to guarantee a unique solution of the eq-system. auto backwardTransitions = epochModel.epochMatrix.transpose(true); @@ -547,7 +547,7 @@ namespace storm { template void MultiDimensionalRewardUnfolding::setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat) { STORM_LOG_ASSERT(model.isOfType(storm::models::ModelType::Dtmc), "Trying to set the equation problem format although the model is not deterministic."); - equationSolverProblemFormatForEpochModel = eqSysFormat; + epochModel.equationSolverProblemFormat = eqSysFormat; } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h index e1b96c64c..1c3851b4f 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -6,6 +6,7 @@ #include "storm/storage/SparseMatrix.h" #include "storm/modelchecker/multiobjective/Objective.h" #include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h" +#include "storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h" #include "storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h" #include "storm/modelchecker/prctl/helper/rewardbounded/Dimension.h" #include "storm/models/sparse/Model.h" @@ -27,17 +28,7 @@ namespace storm { typedef typename EpochManager::EpochClass EpochClass; typedef typename std::conditional>::type SolutionType; - - struct EpochModel { - bool epochMatrixChanged; - storm::storage::SparseMatrix epochMatrix; - storm::storage::BitVector stepChoices; - std::vector stepSolutions; - std::vector> objectiveRewards; - std::vector objectiveRewardFilter; - storm::storage::BitVector epochInStates; - }; - + /* * * @param model The (preprocessed) model @@ -52,7 +43,7 @@ namespace storm { Epoch getStartEpoch(); std::vector getEpochComputationOrder(Epoch const& startEpoch); - EpochModel& setCurrentEpoch(Epoch const& epoch); + EpochModel& setCurrentEpoch(Epoch const& epoch); void setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat); @@ -117,13 +108,10 @@ namespace storm { std::vector epochModelToProductChoiceMap; std::shared_ptr const> productStateToEpochModelInStateMap; std::set possibleEpochSteps; - - EpochModel epochModel; + + EpochModel epochModel; boost::optional currentEpoch; - - // In case of DTMCs we have different options for the equation problem format the epoch model will have. - boost::optional equationSolverProblemFormatForEpochModel; - + EpochManager epochManager; std::vector> dimensions; From f99a24acd22bc3b620c5686852845ff85a45ba0d Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 29 Jan 2019 13:09:32 +0100 Subject: [PATCH 09/20] more work on quantiles. --- .../helper/rewardbounded/QuantileHelper.cpp | 83 +++++++++++++++---- .../helper/rewardbounded/QuantileHelper.h | 8 ++ 2 files changed, 76 insertions(+), 15 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 336c1d299..83cfd7a6a 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -6,7 +6,11 @@ #include #include "storm/models/sparse/Mdp.h" +#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" #include "storm/storage/expressions/ExpressionManager.h" +#include "storm/storage/expressions/Expressions.h" +#include "storm/storage/BitVector.h" +#include "storm/utility/vector.h" #include "storm/logic/ProbabilityOperatorFormula.h" #include "storm/logic/BoundedUntilFormula.h" @@ -21,9 +25,9 @@ namespace storm { // Intentionally left empty } - std::shared_ptr replaceBoundsByGreaterEqZero(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::set const& dimensionsToReplace) { + std::shared_ptr replaceBoundsByGreaterEqZero(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& dimensionsToReplace) { auto const& origBoundedUntil = boundedUntilOperator.getSubformula().asBoundedUntilFormula(); - STORM_LOG_ASSERT(*(--dimensionsToReplace.end()) < origBoundedUntil.getDimension(), "Tried to replace the bound of a dimension that is higher than the number of dimensions of the formula."); + STORM_LOG_ASSERT(dimensionsToReplace.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; std::vector> lowerBounds, upperBounds; std::vector timeBoundReferences; @@ -32,7 +36,17 @@ namespace storm { leftSubformulas.push_back(origBoundedUntil.getLeftSubformula(dim).asSharedPointer()); rightSubformulas.push_back(origBoundedUntil.getRightSubformula(dim).asSharedPointer()); timeBoundReferences.push_back(origBoundedUntil.getTimeBoundReference(dim)); - if (dimensionsToReplace.count(dim) == 0) { + if (dimensionsToReplace.get(dim)) { + storm::expressions::Expression zero; + if (origBoundedUntil.hasLowerBound(dim)) { + zero = origBoundedUntil.getLowerBound(dim).getManager().rational(0.0); + } else { + STORM_LOG_THROW(origBoundedUntil.hasUpperBound(dim), storm::exceptions::InvalidOperationException, "The given bounded until formula has no cost-bound for one dimension."); + zero = origBoundedUntil.getUpperBound(dim).getManager().rational(0.0); + } + lowerBounds.push_back(storm::logic::TimeBound(false, zero)); + upperBounds.push_back(boost::none); + } else { if (origBoundedUntil.hasLowerBound()) { lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); } else { @@ -43,16 +57,6 @@ namespace storm { } else { upperBounds.push_back(boost::none); } - } else { - storm::expressions::Expression zero; - if (origBoundedUntil.hasLowerBound(dim)) { - zero = origBoundedUntil.getLowerBound(dim).getManager().rational(0.0); - } else { - STORM_LOG_THROW(origBoundedUntil.hasUpperBound(dim), storm::exceptions::InvalidOperationException, "The given bounded until formula has no cost-bound for one dimension."); - zero = origBoundedUntil.getUpperBound(dim).getManager().rational(0.0); - } - lowerBounds.push_back(storm::logic::TimeBound(false, zero)); - upperBounds.push_back(boost::none); } } auto newBoundedUntil = std::make_shared(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences); @@ -61,6 +65,26 @@ namespace storm { + template + uint64_t QuantileHelper::getDimension() const { + return quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().getDimension(); + } + + template + storm::storage::BitVector QuantileHelper::getDimensionsForVariable(storm::expressions::Variable const& var) { + 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::computeMultiDimensionalQuantile(Environment const& env) { std::vector> result; @@ -78,8 +102,37 @@ namespace storm { template bool QuantileHelper::computeUnboundedValue(Environment const& env) { - auto unboundedFormula = replaceBoundsByGreaterEqZero(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::set()); - return false; + auto unboundedFormula = replaceBoundsByGreaterEqZero(quantileFormula.getSubformula().asProbabilityOperatorFormula(), storm::storage::BitVector(getDimension(), true)); + MultiDimensionalRewardUnfolding rewardUnfolding(model, unboundedFormula); + + // 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); + STORM_LOG_ASSERT(epochOrder.size() == 1, "unexpected epoch order size."); + // initialize data that will be needed for each epoch + std::vector x, b; + std::unique_ptr> minMaxSolver; +/* + ValueType precision = rewardUnfolding.getRequiredEpochModelPrecision(initEpoch, storm::utility::convertNumber(storm::settings::getModule().getPrecision())); + Environment preciseEnv = env; + preciseEnv.solver().minMax().setPrecision(storm::utility::convertNumber(precision)); + + auto& epochModel = rewardUnfolding.setCurrentEpoch(initEpoch); + // If the epoch matrix is empty we do not need to solve a linear equation system + if (epochModel.epochMatrix.getEntryCount() == 0) { + rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialMdpEpochModel(dir, epochModel)); + } else { + rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialMdpEpochModel(preciseEnv, dir, epochModel, x, b, minMaxSolver, lowerBound, upperBound)); + } + + ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch); + std::cout << "Numeric result is " << numericResult; + return unboundedFormula->getBound().isSatisfied(numericResult); + */ } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h index 276b9faed..915193fbc 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -5,6 +5,10 @@ namespace storm { class Environment; + namespace storage { + class BitVector; + } + namespace modelchecker { namespace helper { namespace rewardbounded { @@ -22,6 +26,10 @@ namespace storm { private: bool computeUnboundedValue(Environment const& env); + uint64_t getDimension() const; + storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var); + + ModelType const& model; storm::logic::QuantileFormula const& quantileFormula; }; From 6e8aef2acc9f227dcd0fb5524d617d6ea73e95a8 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Wed, 30 Jan 2019 12:54:12 +0100 Subject: [PATCH 10/20] Checking formulas with >=0 bound. --- src/storm/logic/Bound.h | 2 +- .../prctl/helper/rewardbounded/Dimension.h | 18 +++++++ .../MultiDimensionalRewardUnfolding.cpp | 50 ++++++++++++------- .../helper/rewardbounded/ProductModel.cpp | 8 ++- .../helper/rewardbounded/QuantileHelper.cpp | 12 +---- .../transformer/EndComponentEliminator.h | 2 +- 6 files changed, 59 insertions(+), 33 deletions(-) diff --git a/src/storm/logic/Bound.h b/src/storm/logic/Bound.h index a39edb5dc..0b152ea06 100644 --- a/src/storm/logic/Bound.h +++ b/src/storm/logic/Bound.h @@ -17,7 +17,7 @@ namespace storm { storm::expressions::Expression threshold; template - bool isSatisfied(ValueType const& compareValue) { + bool isSatisfied(ValueType const& compareValue) const { ValueType thresholdAsValueType = storm::utility::convertNumber(threshold.evaluateAsRational()); switch(comparisonType) { case ComparisonType::Greater: diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h index 76127e5e8..256225616 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -12,13 +12,31 @@ namespace storm { template struct Dimension { + /// The formula describing this dimension std::shared_ptr formula; + + /// The index of the associated objective uint64_t objectiveIndex; + + /// 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 not bounded at all (i.e., it has lower bound >= 0) + bool isNotBounded; + + /// True iff the objective is upper bounded, false if it has a lower bound or no bound at all. bool isUpperBounded; + + /// Multiplying an epoch value with this factor yields the reward/cost in the original domain. ValueType scalingFactor; + + /// The dimensions that are not satisfiable whenever the bound of this dimension is violated storm::storage::BitVector dependentDimensions; + + /// The maximal epoch value that needs to be considered for this dimension boost::optional maxValue; + + /// Whether we minimize/maximize the objective for this dimension boost::optional optimizationDirection; }; } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 4bd08a021..c31e64c80 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -86,15 +86,20 @@ namespace storm { memLabel = "_" + memLabel; } dimension.memoryLabel = memLabel; - dimension.isUpperBounded = subformula.hasUpperBound(dim); - // for simplicity we do not allow intervals or unbounded formulas. - // TODO: Quantiles: allow unbounded formulas - STORM_LOG_THROW(subformula.hasLowerBound(dim) != dimension.isUpperBounded, 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."); + // 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."); - if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { + dimension.isUpperBounded = subformula.hasUpperBound(dim); + // 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()))) { + dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 0)); + dimension.scalingFactor = storm::utility::zero(); + dimension.isNotBounded = true; + } 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.isNotBounded = false; } else { STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); std::string const& rewardName = subformula.getTimeBoundReference(dim).getRewardName(); @@ -105,6 +110,7 @@ namespace storm { auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(actionRewards); dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first)); dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second); + dimension.isNotBounded = false; } dimensions.emplace_back(std::move(dimension)); } @@ -115,6 +121,7 @@ namespace storm { dimension.formula = subformula.restrictToDimension(dim); dimension.objectiveIndex = objIndex; dimension.isUpperBounded = true; + dimension.isNotBounded = false; if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); dimension.scalingFactor = storm::utility::one(); @@ -218,20 +225,23 @@ namespace storm { isStrict = dimFormula.asCumulativeRewardFormula().isBoundStrict(); } - if (bound.containsVariables()) { + if (!bound.containsVariables()) { ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); - STORM_LOG_THROW(dimensions[dim].isUpperBounded || isStrict || !storm::utility::isZero(discretizedBound), storm::exceptions::NotSupportedException, "Lower bounds need to be either strict or greater than zero."); - discretizedBound /= dimensions[dim].scalingFactor; - if (storm::utility::isInteger(discretizedBound)) { - if (isStrict == dimensions[dim].isUpperBounded) { - discretizedBound -= storm::utility::one(); + // 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 needs extra treatment + if (!dimensions[dim].isNotBounded) { + discretizedBound /= dimensions[dim].scalingFactor; + if (storm::utility::isInteger(discretizedBound)) { + if (isStrict == dimensions[dim].isUpperBounded) { + discretizedBound -= storm::utility::one(); + } + } else { + discretizedBound = storm::utility::floor(discretizedBound); } - } else { - discretizedBound = storm::utility::floor(discretizedBound); + uint64_t dimensionValue = storm::utility::convertNumber(discretizedBound); + STORM_LOG_THROW(epochManager.isValidDimensionValue(dimensionValue), storm::exceptions::NotSupportedException, "The bound " << bound << " is too high for the considered number of dimensions."); + dimensions[dim].maxValue = dimensionValue; } - uint64_t dimensionValue = storm::utility::convertNumber(discretizedBound); - STORM_LOG_THROW(epochManager.isValidDimensionValue(dimensionValue), storm::exceptions::NotSupportedException, "The bound " << bound << " is too high for the considered number of dimensions."); - dimensions[dim].maxValue = dimensionValue; } } } @@ -240,8 +250,12 @@ namespace storm { typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch() { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); - epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + if (dimensions[dim].isNotBounded) { + epochManager.setBottomDimension(startEpoch, dim); + } else { + STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); + epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } } STORM_LOG_TRACE("Start epoch is " << epochManager.toString(startEpoch)); return startEpoch; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index af41c8a4f..bc886afcd 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -460,8 +460,12 @@ namespace storm { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); - epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + if (dimensions[dim].isNotBounded) { + epochManager.setBottomDimension(startEpoch, dim); + } else { + STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); + epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } } std::set seenEpochs({startEpoch}); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 83cfd7a6a..5be068cf3 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -116,23 +116,13 @@ namespace storm { // initialize data that will be needed for each epoch std::vector x, b; std::unique_ptr> minMaxSolver; -/* - ValueType precision = rewardUnfolding.getRequiredEpochModelPrecision(initEpoch, storm::utility::convertNumber(storm::settings::getModule().getPrecision())); - Environment preciseEnv = env; - preciseEnv.solver().minMax().setPrecision(storm::utility::convertNumber(precision)); auto& epochModel = rewardUnfolding.setCurrentEpoch(initEpoch); - // If the epoch matrix is empty we do not need to solve a linear equation system - if (epochModel.epochMatrix.getEntryCount() == 0) { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialMdpEpochModel(dir, epochModel)); - } else { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialMdpEpochModel(preciseEnv, dir, epochModel, x, b, minMaxSolver, lowerBound, upperBound)); - } + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch); std::cout << "Numeric result is " << numericResult; return unboundedFormula->getBound().isSatisfied(numericResult); - */ } diff --git a/src/storm/transformer/EndComponentEliminator.h b/src/storm/transformer/EndComponentEliminator.h index 378edc86a..d2489eadc 100644 --- a/src/storm/transformer/EndComponentEliminator.h +++ b/src/storm/transformer/EndComponentEliminator.h @@ -42,7 +42,7 @@ namespace storm { * If addSelfLoopAtSinkStates is true, such rows get a selfloop (with value 1). Otherwise, the row remains empty. */ static EndComponentEliminatorReturnType transform(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& possibleECRows, storm::storage::BitVector const& addSinkRowStates, bool addSelfLoopAtSinkStates = false) { - STORM_LOG_DEBUG("Invoked EndComponentEliminator on matrix with " << originalMatrix.getRowGroupCount() << " row groups."); + STORM_LOG_DEBUG("Invoked EndComponentEliminator on matrix with " << originalMatrix.getRowGroupCount() << " row groups and " << subsystemStates.getNumberOfSetBits() << " subsystem states."); storm::storage::MaximalEndComponentDecomposition ecs = computeECs(originalMatrix, possibleECRows, subsystemStates); From 4ac23d630ffbc6faff71a34a691596a8b5b08873 Mon Sep 17 00:00:00 2001 From: TimQu Date: Thu, 31 Jan 2019 11:53:05 +0100 Subject: [PATCH 11/20] quantiles: Added support for formulas with trivial bounds (i.e., >=0). --- .../prctl/helper/rewardbounded/Dimension.h | 4 ++-- .../prctl/helper/rewardbounded/EpochModel.cpp | 8 ++++---- .../MultiDimensionalRewardUnfolding.cpp | 18 ++++++++---------- .../helper/rewardbounded/ProductModel.cpp | 14 ++++++++------ src/storm/transformer/EndComponentEliminator.h | 2 +- 5 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h index 256225616..34a3bef29 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -21,8 +21,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 not bounded at all (i.e., it has lower bound >= 0) - bool isNotBounded; + /// 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; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp index 37c39e691..e8b059f32 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp @@ -238,10 +238,10 @@ namespace storm { } } - template class EpochModel; - template class EpochModel; - template class EpochModel; - template class EpochModel; + template struct EpochModel; + template struct EpochModel; + template struct EpochModel; + template struct EpochModel; } } } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index c31e64c80..902128890 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -95,11 +95,11 @@ namespace storm { if ((!subformula.hasLowerBound(dim) && !subformula.hasUpperBound(dim)) || (subformula.hasLowerBound(dim) && !subformula.isLowerBoundStrict(dim) && !subformula.getLowerBound(dim).containsVariables() && storm::utility::isZero(subformula.getLowerBound(dim).evaluateAsRational()))) { dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 0)); dimension.scalingFactor = storm::utility::zero(); - dimension.isNotBounded = true; + 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.isNotBounded = false; + dimension.isBounded = true; } else { STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); std::string const& rewardName = subformula.getTimeBoundReference(dim).getRewardName(); @@ -110,7 +110,7 @@ namespace storm { auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(actionRewards); dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first)); dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second); - dimension.isNotBounded = false; + dimension.isBounded = true; } dimensions.emplace_back(std::move(dimension)); } @@ -121,7 +121,7 @@ namespace storm { dimension.formula = subformula.restrictToDimension(dim); dimension.objectiveIndex = objIndex; dimension.isUpperBounded = true; - dimension.isNotBounded = false; + 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(); @@ -229,7 +229,7 @@ namespace storm { ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); // 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 needs extra treatment - if (!dimensions[dim].isNotBounded) { + if (dimensions[dim].isBounded) { discretizedBound /= dimensions[dim].scalingFactor; if (storm::utility::isInteger(discretizedBound)) { if (isStrict == dimensions[dim].isUpperBounded) { @@ -250,11 +250,11 @@ namespace storm { typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch() { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (dimensions[dim].isNotBounded) { - epochManager.setBottomDimension(startEpoch, dim); - } else { + if (dimensions[dim].isBounded) { STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } else { + epochManager.setBottomDimension(startEpoch, dim); } } STORM_LOG_TRACE("Start epoch is " << epochManager.toString(startEpoch)); @@ -745,7 +745,6 @@ namespace storm { } predecessorEpochs.erase(currentEpoch.get()); successorEpochs.erase(currentEpoch.get()); - STORM_LOG_ASSERT(!predecessorEpochs.empty(), "There are no predecessors for the epoch " << epochManager.toString(currentEpoch.get())); // clean up solutions that are not needed anymore for (auto const& successorEpoch : successorEpochs) { @@ -794,7 +793,6 @@ namespace storm { template typename MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::getInitialStateResult(Epoch const& epoch) { STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "The model has multiple initial states."); - STORM_LOG_ASSERT(!epochManager.hasBottomDimension(epoch), "Tried to get the initial state result in an epoch that still contains at least one bottom dimension."); return getStateSolution(epoch, productModel->getInitialProductState(*model.getInitialStates().begin(), model.getInitialStates())); } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index bc886afcd..5072082ee 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -460,11 +460,11 @@ namespace storm { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (dimensions[dim].isNotBounded) { - epochManager.setBottomDimension(startEpoch, dim); - } else { + if (dimensions[dim].isBounded) { STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } else { + epochManager.setBottomDimension(startEpoch, dim); } } @@ -491,22 +491,24 @@ namespace storm { void ProductModel::computeReachableStates(EpochClass const& epochClass, std::vector const& predecessors) { storm::storage::BitVector bottomDimensions(epochManager.getDimensionCount(), false); + bool considerInitialStates = true; for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { if (epochManager.isBottomDimensionEpochClass(epochClass, dim)) { bottomDimensions.set(dim, true); + if (dimensions[dim].isBounded) { + considerInitialStates = false; + } } } storm::storage::BitVector nonBottomDimensions = ~bottomDimensions; storm::storage::BitVector ecInStates(getProduct().getNumberOfStates(), false); - - if (!epochManager.hasBottomDimensionEpochClass(epochClass)) { + if (considerInitialStates) { for (auto const& initState : getProduct().getInitialStates()) { uint64_t transformedInitState = transformProductState(initState, epochClass, memoryStateManager.getInitialMemoryState()); ecInStates.set(transformedInitState, true); } } - for (auto const& predecessor : predecessors) { storm::storage::BitVector positiveStepDimensions(epochManager.getDimensionCount(), false); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { diff --git a/src/storm/transformer/EndComponentEliminator.h b/src/storm/transformer/EndComponentEliminator.h index d2489eadc..c522e1987 100644 --- a/src/storm/transformer/EndComponentEliminator.h +++ b/src/storm/transformer/EndComponentEliminator.h @@ -53,7 +53,7 @@ namespace storm { keptStates.set(stateActionsPair.first, false); } } - STORM_LOG_DEBUG("Found " << ecs.size() << " end components to eliminate. Keeping " << keptStates.getNumberOfSetBits() << " of " << keptStates.size() << " original states."); + STORM_LOG_DEBUG("Found " << ecs.size() << " end components to eliminate. Keeping " << keptStates.getNumberOfSetBits() << " of " << keptStates.size() << " original states plus " << ecs.size() << "new end component states."); EndComponentEliminatorReturnType result; std::vector newRowGroupIndices; From 1d5f2410b57bd71e7a957aecf28a63a8d02f4025 Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 1 Feb 2019 19:04:32 +0100 Subject: [PATCH 12/20] rewardBounded/RewardUnfolding: Allowed the case that not all dimensions have a bound a priori. --- .../helper/rewardbounded/EpochManager.cpp | 9 ++++++ .../prctl/helper/rewardbounded/EpochManager.h | 3 +- .../MultiDimensionalRewardUnfolding.cpp | 29 ++++++++++--------- .../MultiDimensionalRewardUnfolding.h | 11 +++++-- .../helper/rewardbounded/ProductModel.cpp | 20 +++++++++++-- 5 files changed, 53 insertions(+), 19 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp index 8de428483..e248a74da 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp @@ -186,6 +186,15 @@ namespace storm { epoch |= (value << (dimension * bitsPerDimension)); } + void EpochManager::setDimensionOfEpochClass(EpochClass& epochClass, uint64_t const& dimension, bool const& setToBottom) const { + STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count."); + if (setToBottom) { + epochClass |= (1 << dimension); + } else { + epochClass &= ~(1 << dimension); + } + } + bool EpochManager::isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const { STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count."); return (epoch | (dimensionBitMask << (dimension * bitsPerDimension))) == epoch; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h index e4ddffb24..6784cdac4 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h @@ -13,7 +13,7 @@ namespace storm { public: typedef uint64_t Epoch; // The number of reward steps that are "left" for each dimension - typedef uint64_t EpochClass; // The number of reward steps that are "left" for each dimension + typedef uint64_t EpochClass; // Encodes the dimensions of an epoch that are bottom. Two epoch models within the same class have the same graph structure. EpochManager(); EpochManager(uint64_t dimensionCount); @@ -40,6 +40,7 @@ namespace storm { void setBottomDimension(Epoch& epoch, uint64_t const& dimension) const; void setDimensionOfEpoch(Epoch& epoch, uint64_t const& dimension, uint64_t const& value) const; // assumes that the value is valid, i.e., small enough + void setDimensionOfEpochClass(EpochClass& epochClass, uint64_t const& dimension, bool const& setToBottom) const; bool isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const; bool isBottomDimensionEpochClass(EpochClass const& epochClass, uint64_t const& dimension) const; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 902128890..c710d1a04 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -31,7 +31,7 @@ namespace storm { } template - MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula) : model(model) { + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::function&, EpochManager const&)> const& epochStepsCallback) : model(model) { if (objectiveFormula->isProbabilityOperatorFormula()) { if (objectiveFormula->getSubformula().isMultiObjectiveFormula()) { @@ -53,18 +53,28 @@ namespace storm { objective.considersComplementaryEvent = false; objectives.push_back(std::move(objective)); - initialize(); + initialize(epochStepsCallback); } template - void MultiDimensionalRewardUnfolding::initialize() { + void MultiDimensionalRewardUnfolding::initialize(std::function&, EpochManager const&)> const& epochStepsCallback) { maxSolutionsStored = 0; STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives."); std::vector epochSteps; initializeObjectives(epochSteps); + + if (epochStepsCallback) { + epochStepsCallback(epochSteps, epochManager); + } + // collect which epoch steps are possible + possibleEpochSteps.clear(); + for (auto const& step : epochSteps) { + possibleEpochSteps.insert(step); + } + initializeMemoryProduct(epochSteps); } @@ -186,12 +196,6 @@ namespace storm { epochSteps.push_back(step); } - // collect which epoch steps are possible - possibleEpochSteps.clear(); - for (auto const& step : epochSteps) { - possibleEpochSteps.insert(step); - } - // Set the maximal values we need to consider for each dimension computeMaxDimensionValues(); @@ -247,13 +251,13 @@ namespace storm { } template - typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch() { + typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch(bool setUnknownDimsToBottom) { Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (dimensions[dim].isBounded) { - STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); + if (dimensions[dim].isBounded && 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."); epochManager.setBottomDimension(startEpoch, dim); } } @@ -799,7 +803,6 @@ namespace storm { template typename MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) { STORM_LOG_ASSERT(model.getInitialStates().get(initialStateIndex), "The given model state is not an initial state."); - STORM_LOG_ASSERT(!epochManager.hasBottomDimension(epoch), "Tried to get the initial state result in an epoch that still contains at least one bottom dimension."); return getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates())); } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h index 1c3851b4f..77f1ec64b 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -36,11 +36,16 @@ namespace storm { * */ MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::vector> const& objectives); - MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula); + MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula, std::function&, EpochManager const&)> const& epochStepsCallback = nullptr); ~MultiDimensionalRewardUnfolding() = default; - Epoch getStartEpoch(); + /*! + * Retrieves the desired epoch that needs to be analyzed to compute the reward bounded values. + * @param setUnknownDimsToBottom if true, dimensions for which no maxValue is known are set to the bottom dimension. If false, an exception is thrown if there are unknown maxValues. + */ + Epoch getStartEpoch(bool setUnknownDimsToBottom = false); + std::vector getEpochComputationOrder(Epoch const& startEpoch); EpochModel& setCurrentEpoch(Epoch const& epoch); @@ -68,7 +73,7 @@ namespace storm { private: void setCurrentEpochClass(Epoch const& epoch); - void initialize(); + void initialize(std::function&, EpochManager const&)> const& epochStepsCallback = nullptr); void initializeObjectives(std::vector& epochSteps); void computeMaxDimensionValues(); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index 5072082ee..a2c12698b 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -458,10 +458,11 @@ namespace storm { template void ProductModel::collectReachableEpochClasses(std::set>& reachableEpochClasses, std::set const& possibleSteps) const { + // Get the start epoch according to the given bounds. + // For dimensions for which no bound (aka. maxValue) is known, we will later overapproximate the set of reachable classes. Epoch startEpoch = epochManager.getZeroEpoch(); for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { - if (dimensions[dim].isBounded) { - STORM_LOG_ASSERT(dimensions[dim].maxValue, "No max-value for dimension " << dim << " was given."); + if (dimensions[dim].maxValue) { epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); } else { epochManager.setBottomDimension(startEpoch, dim); @@ -485,6 +486,21 @@ namespace storm { } } } + + // Also treat dimensions without a priori bound + for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { + if (!dimensions[dim].maxValue) { + std::vector newClasses; + for (auto const& c : reachableEpochClasses) { + auto newClass = c; + epochManager.setDimensionOfEpochClass(newClass, dim, false); + newClasses.push_back(newClass); + } + for (auto const& c: newClasses) { + reachableEpochClasses.insert(c); + } + } + } } template From a99dd5e4d1d73e9ca87757c5c1e7d29e9a48944d Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 1 Feb 2019 19:05:57 +0100 Subject: [PATCH 13/20] quantiles: Better code re-usage, better structure, support for 'open' and 'non-open' dimensions, single dimensional quantiles should work now. --- .../prctl/SparseMdpPrctlModelChecker.cpp | 2 +- .../helper/rewardbounded/QuantileHelper.cpp | 238 ++++++++++++++++-- .../helper/rewardbounded/QuantileHelper.h | 31 ++- 3 files changed, 236 insertions(+), 35 deletions(-) diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index b16bd7b3c..759d7a976 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -233,7 +233,7 @@ namespace storm { uint64_t initialState = *this->getModel().getInitialStates().begin(); helper::rewardbounded::QuantileHelper qHelper(this->getModel(), checkTask.getFormula()); - auto res = qHelper.computeMultiDimensionalQuantile(env); + auto res = qHelper.computeQuantile(env); if (res.size() == 1 && res.front().size() == 1) { return std::unique_ptr(new ExplicitQuantitativeCheckResult(initialState, std::move(res.front().front()))); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 5be068cf3..d20f4a7da 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -4,12 +4,14 @@ #include #include #include +#include #include "storm/models/sparse/Mdp.h" #include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" #include "storm/storage/expressions/ExpressionManager.h" #include "storm/storage/expressions/Expressions.h" #include "storm/storage/BitVector.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" #include "storm/utility/vector.h" #include "storm/logic/ProbabilityOperatorFormula.h" @@ -23,11 +25,23 @@ namespace storm { template QuantileHelper::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) { // Intentionally left empty + /* Todo: Current assumption: + * Subformula is always prob operator with bounded until + * Each bound variable occurs at most once (change this?) + * cost bounds can either be <= or > (the epochs returned by the reward unfolding assume this. One could translate them, though) + * 'reasonable' quantile (e.g. not quantile(max B, Pmax>0.5 [F <=B G]) (we could also filter out these things early on) + */ } - std::shared_ptr replaceBoundsByGreaterEqZero(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& dimensionsToReplace) { + enum class BoundTransformation { + None, + GreaterZero, + GreaterEqualZero, + LessEqualZero + }; + std::shared_ptr transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector const& transformations) { auto const& origBoundedUntil = boundedUntilOperator.getSubformula().asBoundedUntilFormula(); - STORM_LOG_ASSERT(dimensionsToReplace.size() == origBoundedUntil.getDimension(), "Tried to replace the bound of a dimension that is higher than the number of dimensions of the formula."); + 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; std::vector> lowerBounds, upperBounds; std::vector timeBoundReferences; @@ -36,7 +50,19 @@ namespace storm { leftSubformulas.push_back(origBoundedUntil.getLeftSubformula(dim).asSharedPointer()); rightSubformulas.push_back(origBoundedUntil.getRightSubformula(dim).asSharedPointer()); timeBoundReferences.push_back(origBoundedUntil.getTimeBoundReference(dim)); - if (dimensionsToReplace.get(dim)) { + if (transformations[dim] == BoundTransformation::None) { + if (origBoundedUntil.hasLowerBound()) { + lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); + } else { + lowerBounds.push_back(boost::none); + } + if (origBoundedUntil.hasUpperBound()) { + upperBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim))); + } else { + upperBounds.push_back(boost::none); + } + } else { + // We need a zero expression in all other cases storm::expressions::Expression zero; if (origBoundedUntil.hasLowerBound(dim)) { zero = origBoundedUntil.getLowerBound(dim).getManager().rational(0.0); @@ -44,17 +70,12 @@ namespace storm { STORM_LOG_THROW(origBoundedUntil.hasUpperBound(dim), storm::exceptions::InvalidOperationException, "The given bounded until formula has no cost-bound for one dimension."); zero = origBoundedUntil.getUpperBound(dim).getManager().rational(0.0); } - lowerBounds.push_back(storm::logic::TimeBound(false, zero)); - upperBounds.push_back(boost::none); - } else { - if (origBoundedUntil.hasLowerBound()) { - lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); - } else { + if (transformations[dim] == BoundTransformation::LessEqualZero) { lowerBounds.push_back(boost::none); - } - if (origBoundedUntil.hasUpperBound()) { - lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isUpperBoundStrict(dim), origBoundedUntil.getUpperBound(dim))); + upperBounds.push_back(storm::logic::TimeBound(false, zero)); } else { + STORM_LOG_ASSERT(transformations[dim] == BoundTransformation::GreaterZero || transformations[dim] == BoundTransformation::GreaterEqualZero, "Unhandled bound transformation."); + lowerBounds.push_back(storm::logic::TimeBound(transformations[dim] == BoundTransformation::GreaterZero, zero)); upperBounds.push_back(boost::none); } } @@ -63,15 +84,36 @@ namespace storm { return std::make_shared(newBoundedUntil, boundedUntilOperator.getOperatorInformation()); } - - template uint64_t QuantileHelper::getDimension() const { return quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().getDimension(); } template - storm::storage::BitVector QuantileHelper::getDimensionsForVariable(storm::expressions::Variable const& var) { + storm::storage::BitVector QuantileHelper::getOpenDimensions() const { + auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); + storm::storage::BitVector res(getDimension()); + for (uint64_t dim = 0; dim < getDimension(); ++dim) { + res.set(dim, boundedUntil.hasLowerBound(dim) ? boundedUntil.getLowerBound(dim).containsVariables() : boundedUntil.getUpperBound(dim).containsVariables()); + } + 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(); + for (auto const& boundVar : quantileFormula.getBoundVariables()) { + if (boundVar.second == dimVar) { + return boundVar.first; + } + } + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "The bound variable '" << dimVar.getName() << "' is not specified within the quantile formula '" << quantileFormula << "'."); + return quantileFormula.getOptimizationDirection(); + } + + template + 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) { @@ -86,24 +128,166 @@ namespace storm { } template - std::vector> QuantileHelper::computeMultiDimensionalQuantile(Environment const& env) { + std::vector> QuantileHelper::computeQuantile(Environment const& env) { std::vector> result; - if (!computeUnboundedValue(env)) { + if (getOpenDimensions().getNumberOfSetBits() == 1) { + 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) { + auto quantileRes = computeQuantileForDimension(env, dimension); + result = {{storm::utility::convertNumber(quantileRes.first) * quantileRes.second}}; + } else if (maxProbSatisfiesFormula) { + // i.e., all bound values satisfy the formula + bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension)); + bool upperCostBound = boundedUntilOperator.getSubformula().asBoundedUntilFormula().hasUpperBound(dimension); + if (minimizingRewardBound) { + result = {{ (upperCostBound ? storm::utility::zero() : -storm::utility::one())}}; + } else { + result = {{ storm::utility::infinity()}}; + } + } else { + // i.e., no bound value satisfies the formula + result = {{}}; + } + } else if (getOpenDimensions().getNumberOfSetBits() == 2) { + + } 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 result; + /* + return unboundedFormula->getBound().isSatisfied(numericResult); + + if (!checkUnboundedValue(env)) { STORM_LOG_INFO("Bound is not satisfiable."); result.emplace_back(); for (auto const& bv : quantileFormula.getBoundVariables()) { result.front().push_back(storm::solver::minimize(bv.first) ? storm::utility::infinity() : -storm::utility::infinity()); } } else { + std::vector limitProbCheckResults; + for (uint64_t dim = 0; dim < getDimension(); ++dim) { + storm::storage::BitVector dimAsBitVector(getDimension(), false); + dimAsBitVector.set(dim, true); + limitProbCheckResults.push_back(checkLimitProbability(env, dimAsBitVector)); + auto quantileRes = computeQuantileForDimension(env, dim); + std::cout << "Quantile for dim " << dim << " is " << (storm::utility::convertNumber(quantileRes.first) * quantileRes.second) << std::endl; + } + result = {{27}}; } - return result; + return result;*/ } template - bool QuantileHelper::computeUnboundedValue(Environment const& env) { - auto unboundedFormula = replaceBoundsByGreaterEqZero(quantileFormula.getSubformula().asProbabilityOperatorFormula(), storm::storage::BitVector(getDimension(), true)); - MultiDimensionalRewardUnfolding rewardUnfolding(model, unboundedFormula); + std::pair QuantileHelper::computeQuantileForDimension(Environment const& env, uint64_t dimension) const { + // We assume that their is one bound value that violates the quantile and one bound value that satisfies it. + + // 'Drop' all other open bounds + std::vector bts(getDimension(), BoundTransformation::None); + for (auto const& d : getOpenDimensions()) { + if (d != dimension) { + bts[d] = BoundTransformation::GreaterEqualZero; + } + } + auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts); + MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula); + + bool upperCostBound = transformedFormula->getSubformula().asBoundedUntilFormula().hasUpperBound(dimension); + // bool lowerProbabilityBound = storm::logic::isLowerBound(transformedFormula->getBound().comparisonType); + 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; + + int64_t epochValue = 0; // -1 is actally allowed since we express >=0 as >-1 + 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)); + exploredEpochs.insert(epoch); + } + } + 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); + // 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) { + break; + } else if (!minimizingRewardBound && !propertySatisfied) { + STORM_LOG_ASSERT(epochValue > 0 || !upperCostBound, "The property does not seem to be satisfiable. This case should have been treated earlier."); + --epochValue; + break; + } + ++epochValue; + } + return {epochValue, rewardUnfolding.getDimension(dimension).scalingFactor}; + } + + 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 + + std::vector bts(getDimension(), BoundTransformation::None); + storm::storage::BitVector minimizingLowerBoundedDimensions(getDimension(), false); + + for (auto const& d : getOpenDimensions()) { + if (minimizingDimensions.get(d)) { + bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d); + 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."); + } + } 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); + } + } + } + } + + 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); + } + } + } + }); // Get lower and upper bounds for the solution. auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); @@ -112,20 +296,20 @@ namespace storm { // Initialize epoch models auto initEpoch = rewardUnfolding.getStartEpoch(); auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch); - STORM_LOG_ASSERT(epochOrder.size() == 1, "unexpected epoch order size."); // initialize data that will be needed for each epoch std::vector x, b; std::unique_ptr> minMaxSolver; - auto& epochModel = rewardUnfolding.setCurrentEpoch(initEpoch); - rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); + for (auto const& epoch : epochOrder) { + auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); + } ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch); - std::cout << "Numeric result is " << numericResult; - return unboundedFormula->getBound().isSatisfied(numericResult); + STORM_LOG_TRACE("Extremal probability for minimizing dimensions " << minimizingDimensions << " is " << numericResult << "."); + return transformedFormula->getBound().isSatisfied(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 915193fbc..789e1d727 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -19,16 +19,33 @@ namespace storm { public: QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula); - std::vector> computeMultiDimensionalQuantile(Environment const& env); - - + std::vector> computeQuantile(Environment const& env); private: - bool computeUnboundedValue(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. + */ + ValueType computeExtremalValue(Environment const& env, storm::storage::BitVector const& minimizingDimensions) const; + + + /// Computes the quantile with respect to the given dimension + std::pair computeQuantileForDimension(Environment const& env, uint64_t dim) const; + + /*! + * Gets the number of dimensions of the underlying boudned until formula + */ uint64_t getDimension() const; - storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var); - + + /*! + * Gets the dimensions that are open, i.e., for which the bound value is not fixed + * @return + */ + storm::storage::BitVector getOpenDimensions() const; + + storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var) const; + storm::solver::OptimizationDirection const& getOptimizationDirForDimension(uint64_t const& dim) const; ModelType const& model; storm::logic::QuantileFormula const& quantileFormula; From 6aeb75e3bddfcbe73ced9bdbac47cdadfcedb45a Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 5 Feb 2019 19:14:15 +0100 Subject: [PATCH 14/20] quantiles: Supporting two-dimensional quantiles with the same optimization direction of quantile bounds (max,max or min,min). --- .../MultiDimensionalRewardUnfolding.cpp | 4 +- .../helper/rewardbounded/ProductModel.cpp | 2 +- .../helper/rewardbounded/QuantileHelper.cpp | 288 +++++++++++++++--- .../helper/rewardbounded/QuantileHelper.h | 1 + 4 files changed, 255 insertions(+), 40 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index c710d1a04..1d57dd2b7 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -230,10 +230,10 @@ namespace storm { } if (!bound.containsVariables()) { - ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); // 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 needs extra treatment + // Thus, >=N would become >N-1. However, note that the case N=0 is treated separately. if (dimensions[dim].isBounded) { + ValueType discretizedBound = storm::utility::convertNumber(bound.evaluateAsRational()); discretizedBound /= dimensions[dim].scalingFactor; if (storm::utility::isInteger(discretizedBound)) { if (isStrict == dimensions[dim].isUpperBounded) { diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index a2c12698b..57fbf4980 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -511,7 +511,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) { + if (dimensions[dim].isBounded && dimensions[dim].maxValue) { considerInitialStates = false; } } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index d20f4a7da..a06907e2f 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include "storm/models/sparse/Mdp.h" #include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" @@ -17,6 +16,8 @@ #include "storm/logic/ProbabilityOperatorFormula.h" #include "storm/logic/BoundedUntilFormula.h" +#include "storm/exceptions/NotSupportedException.h" + namespace storm { namespace modelchecker { namespace helper { @@ -25,10 +26,20 @@ namespace storm { template QuantileHelper::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) { // Intentionally left empty - /* Todo: Current assumption: + 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.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(); + + // 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 + + /* Todo: Current assumptions: * Subformula is always prob operator with bounded until * Each bound variable occurs at most once (change this?) - * cost bounds can either be <= or > (the epochs returned by the reward unfolding assume this. One could translate them, though) + * 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) */ } @@ -141,10 +152,8 @@ namespace storm { result = {{storm::utility::convertNumber(quantileRes.first) * quantileRes.second}}; } else if (maxProbSatisfiesFormula) { // i.e., all bound values satisfy the formula - bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension)); - bool upperCostBound = boundedUntilOperator.getSubformula().asBoundedUntilFormula().hasUpperBound(dimension); - if (minimizingRewardBound) { - result = {{ (upperCostBound ? storm::utility::zero() : -storm::utility::one())}}; + if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) { + result = {{ storm::utility::zero() }}; } else { result = {{ storm::utility::infinity()}}; } @@ -153,52 +162,246 @@ namespace storm { result = {{}}; } } else if (getOpenDimensions().getNumberOfSetBits() == 2) { - + result = computeTwoDimensionalQuantile(env); } 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 result; - /* - return unboundedFormula->getBound().isSatisfied(numericResult); - - if (!checkUnboundedValue(env)) { - STORM_LOG_INFO("Bound is not satisfiable."); - result.emplace_back(); - for (auto const& bv : quantileFormula.getBoundVariables()) { - result.front().push_back(storm::solver::minimize(bv.first) ? storm::utility::infinity() : -storm::utility::infinity()); + } + + 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; + } } - } else { - std::vector limitProbCheckResults; - for (uint64_t dim = 0; dim < getDimension(); ++dim) { - storm::storage::BitVector dimAsBitVector(getDimension(), false); - dimAsBitVector.set(dim, true); - limitProbCheckResults.push_back(checkLimitProbability(env, dimAsBitVector)); - auto quantileRes = computeQuantileForDimension(env, dim); - std::cout << "Quantile for dim " << dim << " is " << (storm::utility::convertNumber(quantileRes.first) * quantileRes.second) << std::endl; + if (!p1Dominated) { + result.push_back(p1); } - - result = {{27}}; } - return result;*/ + points = std::move(result); } + + template + std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment const& env) { + std::vector> result; + auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); + storm::logic::ComparisonType probabilityBound = probOpFormula.getBound().comparisonType; + 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 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) { + 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])); + 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 + 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); + // 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::cout << "Found starting point. Beginning 2D exploration" << std::endl; + + 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; + + std::vector 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 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 convertedP; + for (uint64_t i = 0; i < 2; ++i) { + if (p[i] == infinity) { + convertedP.push_back(storm::utility::infinity()); + } 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(p[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor); + } + } + if (!convertedP.empty()) { + result.push_back(convertedP); + } + } + filterDominatedPoints(result, optimizationDirections); + } else if (maxmaxProbSatisfiesFormula) { + // 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 std::pair QuantileHelper::computeQuantileForDimension(Environment const& env, uint64_t dimension) const { - // We assume that their is one bound value that violates the quantile and one bound value that satisfies it. + // We assume that there is one bound value that violates the quantile and one bound value that satisfies it. - // 'Drop' all other open bounds + // Let all other open bounds approach infinity std::vector bts(getDimension(), BoundTransformation::None); + storm::storage::BitVector otherLowerBoundedDimensions(getDimension(), false); for (auto const& d : getOpenDimensions()) { if (d != dimension) { - bts[d] = BoundTransformation::GreaterEqualZero; + if (quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasLowerBound(d)) { + bts[d] = BoundTransformation::GreaterZero; + otherLowerBoundedDimensions.set(d, true); + } else { + bts[d] = BoundTransformation::GreaterEqualZero; + } } } + + 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); + 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); + } + } + } + }); - bool upperCostBound = transformedFormula->getSubformula().asBoundedUntilFormula().hasUpperBound(dimension); - // bool lowerProbabilityBound = storm::logic::isLowerBound(transformedFormula->getBound().comparisonType); - bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension)); // initialize data that will be needed for each epoch auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); @@ -206,7 +409,7 @@ namespace storm { std::vector x, b; std::unique_ptr> minMaxSolver; - int64_t epochValue = 0; // -1 is actally allowed since we express >=0 as >-1 + uint64_t epochValue = 0; std::set exploredEpochs; while (true) { auto currEpoch = rewardUnfolding.getStartEpoch(true); @@ -225,10 +428,20 @@ namespace storm { // 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) { - STORM_LOG_ASSERT(epochValue > 0 || !upperCostBound, "The property does not seem to be satisfiable. This case should have been treated earlier."); - --epochValue; + // 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; @@ -278,6 +491,7 @@ namespace storm { } } } + std::cout << "nonmec choices are " << nonMecChoices << std::endl; MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, [&](std::vector& epochSteps, EpochManager const& epochManager) { if (!minimizingLowerBoundedDimensions.empty()) { diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h index 789e1d727..a0cd22fdb 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -22,6 +22,7 @@ namespace storm { std::vector> computeQuantile(Environment const& env); private: + std::vector> computeTwoDimensionalQuantile(Environment const& env); /*! * Computes the infimum over bound values in minimizing dimensions / the supremum over bound values in maximizing dimensions From dd93b1dae9ed84cfff63e277af24b58699e20300 Mon Sep 17 00:00:00 2001 From: TimQu Date: Thu, 7 Feb 2019 16:30:43 +0100 Subject: [PATCH 15/20] rewardbounded: Improved code structure. --- ...ewardBoundedMdpPcaaWeightVectorChecker.cpp | 2 +- .../prctl/helper/SparseDtmcPrctlHelper.cpp | 2 +- .../prctl/helper/SparseMdpPrctlHelper.cpp | 2 +- .../prctl/helper/rewardbounded/Dimension.h | 15 +- .../MultiDimensionalRewardUnfolding.cpp | 168 ++++++++++++----- .../MultiDimensionalRewardUnfolding.h | 21 ++- .../helper/rewardbounded/ProductModel.cpp | 58 ++++-- .../prctl/helper/rewardbounded/ProductModel.h | 6 +- .../helper/rewardbounded/QuantileHelper.cpp | 175 +++++++----------- .../helper/rewardbounded/QuantileHelper.h | 6 +- 10 files changed, 270 insertions(+), 185 deletions(-) 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; From fb7078770dee6304b7d506201c9ec7e1050b2b6b Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 8 Feb 2019 15:32:22 +0100 Subject: [PATCH 16/20] rewardbounded: Various fixes. --- .../MultiDimensionalRewardUnfolding.cpp | 26 +- .../MultiDimensionalRewardUnfolding.h | 12 +- .../helper/rewardbounded/ProductModel.cpp | 8 +- .../prctl/helper/rewardbounded/ProductModel.h | 2 +- .../helper/rewardbounded/QuantileHelper.cpp | 317 ++++++++++++------ .../helper/rewardbounded/QuantileHelper.h | 19 +- 6 files changed, 247 insertions(+), 137 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 450b04b8c..4d7d1deb4 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/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 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())) { @@ -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 - 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> @@ -865,10 +853,7 @@ 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(); - } + 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 + 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 6513b1939..510cdcfd5 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -77,6 +77,12 @@ namespace storm { 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: @@ -89,12 +95,6 @@ namespace storm { 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 4ed49fa1a..aa387d084 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -92,9 +92,8 @@ namespace storm { computeReachableStatesInEpochClasses(); } - template - storm::storage::MemoryStructure ProductModel::computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives) const { + storm::storage::MemoryStructure ProductModel::computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives) { storm::modelchecker::SparsePropositionalModelChecker> mc(model); @@ -649,6 +648,11 @@ namespace storm { return getProductState(getModelState(productState), transformMemoryState(getMemoryState(productState), epochClass, predecessorMemoryState)); } + template + storm::storage::BitVector const& ProductModel::getProb1Objectives() { + return prob1Objectives; + } + template class ProductModel; template class ProductModel; } diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h index 92036761e..6ae408930 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h @@ -52,7 +52,7 @@ namespace storm { private: - storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives) const; + storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives); std::vector computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 84a28a68c..1e6897ee5 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -5,6 +5,8 @@ #include #include +#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(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("0.1"); + env.solver().setLinearEquationSolverPrecision(static_cast(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 + std::pair getLowerUpperBound(storm::Environment const& env, ValueType const& value, bool minMax = true) { + ValueType prec; + bool relative; + if (minMax) { + prec = storm::utility::convertNumber(env.solver().minMax().getPrecision()); + relative = env.solver().minMax().getRelativeTerminationCriterion(); + } else { + prec = storm::utility::convertNumber(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get()); + relative = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).second.get(); + } + if (relative) { + return std::make_pair(value * (1/(prec + 1)), value * (1 + prec/(prec +1))); + } else { + return std::make_pair(value - prec, value + prec); + } + } + template uint64_t QuantileHelper::getDimension() const { return quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().getDimension(); @@ -158,15 +187,21 @@ namespace storm { template std::vector> QuantileHelper::computeQuantile(Environment const& env) { std::vector> 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(quantileRes.first) * quantileRes.second}}; + while (true) { + auto quantileRes = computeQuantileForDimension(envCpy, dimension); + if (quantileRes) { + result = {{storm::utility::convertNumber(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 - std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment const& env) { + std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment& env) { std::vector> 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> 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 - 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 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; - - std::vector 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 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 convertedP; - for (uint64_t i = 0; i < 2; ++i) { - if (p[i] == infinity) { - convertedP.push_back(storm::utility::infinity()); - } 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(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 newResultPoint(2); - newResultPoint[i] = storm::utility::convertNumber(singleQuantileValues[i].first - 1) * singleQuantileValues[i].second; - newResultPoint[1-i] = storm::utility::infinity(); - result.push_back(newResultPoint); - } + 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); } - 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 - std::pair QuantileHelper::computeQuantileForDimension(Environment const& env, uint64_t dimension) const { + 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; + + 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 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 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()) { + 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 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 @@ -391,7 +463,6 @@ namespace storm { 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)); @@ -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 + 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); + } + } else { + return probabilityBound.isSatisfied(numericResult); + } + } } template @@ -463,7 +566,11 @@ namespace storm { 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(); @@ -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>; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h index ce05b5f5a..5cd6bfc8a 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/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> computeQuantile(Environment const& env); private: - std::vector> computeTwoDimensionalQuantile(Environment const& env); - + std::vector> computeTwoDimensionalQuantile(Environment& env); + bool exploreTwoDimensionalQuantile(Environment const& env, std::vector> const& startEpochValues, std::vector& currentEpochValues, std::vector>& 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 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> computeQuantileForDimension(Environment const& env, uint64_t dim) const; /*! * Gets the number of dimensions of the underlying boudned until formula From 5ad100e6526a0f64ea5e1931a0b8147d3d5f5e41 Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 8 Feb 2019 17:00:39 +0100 Subject: [PATCH 17/20] quantiles: Added some statistics. --- .../helper/rewardbounded/QuantileHelper.cpp | 26 +++++++++++-------- .../helper/rewardbounded/QuantileHelper.h | 8 ++++-- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp index 1e6897ee5..3a5ca0239 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -14,6 +14,8 @@ #include "storm/storage/BitVector.h" #include "storm/storage/MaximalEndComponentDecomposition.h" #include "storm/utility/vector.h" +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/CoreSettings.h" #include "storm/logic/ProbabilityOperatorFormula.h" #include "storm/logic/BoundedUntilFormula.h" @@ -185,9 +187,10 @@ namespace storm { } template - std::vector> QuantileHelper::computeQuantile(Environment const& env) { + std::vector> QuantileHelper::computeQuantile(Environment const& env) const { 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(); @@ -201,6 +204,7 @@ namespace storm { break; } increasePrecision(envCpy); + ++numPrecisionRefinements; } } else if (zeroSatisfiesFormula) { // thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula @@ -218,6 +222,10 @@ namespace storm { } else { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The quantile formula considers " << getOpenDimensions().getNumberOfSetBits() << " open dimensions. Only one- or two-dimensional quantiles are supported."); } + if (storm::settings::getModule().isShowStatisticsSet()) { + std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl; + std::cout << "Number of required precision refinements: " << numPrecisionRefinements << std::endl; + } return result; } @@ -253,11 +261,10 @@ namespace storm { } template - std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment& env) { + std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment& env) const { std::vector> result; auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); - storm::logic::ComparisonType probabilityBound = probOpFormula.getBound().comparisonType; auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula(); std::vector dimensionsAsBitVector; std::vector dimensions; @@ -282,15 +289,12 @@ namespace storm { 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]); - 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 is set to infinity - std::cout << "Computing quantile for single dimension " << dimensions[i] << std::endl; 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. // Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i] @@ -302,10 +306,10 @@ namespace storm { --singleQuantileValues[i].first; } } - std::cout << "Found starting point. Beginning 2D exploration" << std::endl; 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 @@ -354,7 +358,6 @@ namespace storm { std::vector x, b; std::unique_ptr> 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]); @@ -379,11 +382,11 @@ namespace storm { 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); - 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); @@ -419,7 +422,6 @@ namespace storm { convertedPoint.push_back(storm::utility::convertNumber(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)); } @@ -482,11 +484,11 @@ namespace storm { 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); - std::cout << "Numeric result for epoch " << rewardUnfolding.getEpochManager().toString(currEpoch) << " is " << currValue << std::endl; bool propertySatisfied; if (env.solver().isForceSoundness()) { @@ -537,6 +539,7 @@ namespace storm { return lowerSat; } else { increasePrecision(env); + ++numPrecisionRefinements; } } else { return probabilityBound.isSatisfied(numericResult); @@ -585,6 +588,7 @@ namespace storm { 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); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h index 5cd6bfc8a..13676a644 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -20,10 +20,10 @@ namespace storm { public: QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula); - std::vector> computeQuantile(Environment const& env); + std::vector> computeQuantile(Environment const& env) const; private: - std::vector> computeTwoDimensionalQuantile(Environment& env); + std::vector> computeTwoDimensionalQuantile(Environment& env) const; bool exploreTwoDimensionalQuantile(Environment const& env, std::vector> const& startEpochValues, std::vector& currentEpochValues, std::vector>& resultPoints) const; /*! @@ -59,6 +59,10 @@ namespace storm { ModelType const& model; storm::logic::QuantileFormula const& quantileFormula; + + /// Statistics + mutable uint64_t numCheckedEpochs; + mutable uint64_t numPrecisionRefinements; }; } } From 04082fb2d62df7b0a239e59d33e04ef24903d651 Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 26 Feb 2019 15:10:25 +0100 Subject: [PATCH 18/20] Added a method to check whether a graph contains a cycle. --- src/storm/utility/graph.cpp | 47 +++++++++++++++++++++++++++++++++++-- src/storm/utility/graph.h | 9 +++++++ 2 files changed, 54 insertions(+), 2 deletions(-) diff --git a/src/storm/utility/graph.cpp b/src/storm/utility/graph.cpp index d17d95edc..fa23afe7f 100644 --- a/src/storm/utility/graph.cpp +++ b/src/storm/utility/graph.cpp @@ -116,6 +116,43 @@ namespace storm { return result; } + template + bool hasCycle(storm::storage::SparseMatrix const& transitionMatrix, boost::optional const& subsystem) { + storm::storage::BitVector unexploredStates; // States that have not been visited yet + storm::storage::BitVector acyclicStates; // States that are known to not lie on a cycle consisting of subsystem states + if (subsystem) { + unexploredStates = subsystem.get(); + acyclicStates = ~subsystem.get(); + } else { + unexploredStates.resize(transitionMatrix.getRowGroupCount(), true); + acyclicStates.resize(transitionMatrix.getRowGroupCount(), false); + } + std::vector dfsStack; + for (uint64_t start = unexploredStates.getNextSetIndex(0); start < unexploredStates.size(); start = unexploredStates.getNextSetIndex(start + 1)) { + dfsStack.push_back(start); + while (!dfsStack.empty()) { + uint64_t state = dfsStack.back(); + if (unexploredStates.get(state)) { + unexploredStates.set(state, false); + for (auto const& entry : transitionMatrix.getRowGroup(start)) { + if (unexploredStates.get(entry.getColumn())) { + dfsStack.push_back(entry.getColumn()); + } else { + if (!acyclicStates.get(entry.getColumn())) { + // The state has been visited before but is not known to be acyclic. + return true; + } + } + } + } else { + acyclicStates.set(state, true); + dfsStack.pop_back(); + } + } + } + return false; + } + template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices) { @@ -1650,6 +1687,8 @@ namespace storm { template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix const& transitionMatrix); + template bool hasCycle(storm::storage::SparseMatrix const& transitionMatrix, boost::optional const& subsystem); + template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices); template std::vector getDistances(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional const& subsystem); @@ -1727,8 +1766,10 @@ namespace storm { template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& constraintStates, storm::storage::BitVector const& targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional const& choiceFilter); template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix const& transitionMatrix); - - template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices); + + template bool hasCycle(storm::storage::SparseMatrix const& transitionMatrix, boost::optional const& subsystem); + + template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices); template std::vector getDistances(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional const& subsystem); @@ -1785,6 +1826,8 @@ namespace storm { template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix const& transitionMatrix); + template bool hasCycle(storm::storage::SparseMatrix const& transitionMatrix, boost::optional const& subsystem); + template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices); template std::vector getDistances(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional const& subsystem); diff --git a/src/storm/utility/graph.h b/src/storm/utility/graph.h index 152cf10b9..6f65bbcfa 100644 --- a/src/storm/utility/graph.h +++ b/src/storm/utility/graph.h @@ -81,6 +81,15 @@ namespace storm { template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix const& transitionMatrix); + + /*! + * Returns true if the graph represented by the given matrix has a cycle + * @param transitionMatrix + * @param subsystem if given, only states in the subsystem are considered for the check. + */ + template + bool hasCycle(storm::storage::SparseMatrix const& transitionMatrix, boost::optional const& subsystem = boost::none); + /*! * Checks whether there is an End Component that * 1. contains at least one of the specified choices and From 88ee0bbf67ed5cc2262afd56adc8739948a74f3e Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 26 Feb 2019 15:11:25 +0100 Subject: [PATCH 19/20] RewardUnfolding: If statistics are enabled, Log when an acyclic epoch model is found. --- .../rewardbounded/MultiDimensionalRewardUnfolding.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index 4d7d1deb4..3b1f36791 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -7,6 +7,8 @@ #include "storm/utility/macros.h" #include "storm/logic/Formulas.h" +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/CoreSettings.h" #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" @@ -377,6 +379,11 @@ namespace storm { if (!currentEpoch || !epochManager.compareEpochClass(epoch, currentEpoch.get())) { setCurrentEpochClass(epoch); epochModel.epochMatrixChanged = true; + if (storm::settings::getModule().isShowStatisticsSet()) { + if (storm::utility::graph::hasCycle(epochModel.epochMatrix)) { + std::cout << "Epoch model for epoch " << epochManager.toString(epoch) << " is cyclic." << std::endl; + } + } } else { epochModel.epochMatrixChanged = false; } @@ -690,7 +697,6 @@ namespace storm { template ValueType MultiDimensionalRewardUnfolding::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) { - // if (storm::settings::getModule().isSoundSet()) { uint64_t sumOfDimensions = 0; for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) { sumOfDimensions += epochManager.getDimensionOfEpoch(startEpoch, dim) + 1; From e89cbf2886ceeaf4c7f9826b9ae7688d87e0b1af Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 27 Feb 2019 00:55:46 +0100 Subject: [PATCH 20/20] fixed cyclic check. --- src/storm/utility/graph.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/storm/utility/graph.cpp b/src/storm/utility/graph.cpp index fa23afe7f..1515a31f7 100644 --- a/src/storm/utility/graph.cpp +++ b/src/storm/utility/graph.cpp @@ -134,7 +134,7 @@ namespace storm { uint64_t state = dfsStack.back(); if (unexploredStates.get(state)) { unexploredStates.set(state, false); - for (auto const& entry : transitionMatrix.getRowGroup(start)) { + for (auto const& entry : transitionMatrix.getRowGroup(state)) { if (unexploredStates.get(entry.getColumn())) { dfsStack.push_back(entry.getColumn()); } else {