diff --git a/CHANGELOG.md b/CHANGELOG.md index 226abe92e..531e57a31 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ The releases of major and minor versions contain an overview of changes since th Version 1.3.x ------------- +### Version 1.3.1 (under development) +- Added support for multi-dimensional quantile queries +- Fixed sparse bisimulation of MDPs (which failed if all non-absorbing states in the quotient are initial) ### Version 1.3.0 (2018/12) - Slightly improved scheduler extraction diff --git a/src/storm-pars/parser/ParameterRegionParser.cpp b/src/storm-pars/parser/ParameterRegionParser.cpp index f6345eef1..c7c5bf98e 100644 --- a/src/storm-pars/parser/ParameterRegionParser.cpp +++ b/src/storm-pars/parser/ParameterRegionParser.cpp @@ -1,3 +1,4 @@ +#include #include "storm-pars/parser/ParameterRegionParser.h" #include "storm/utility/macros.h" @@ -10,13 +11,13 @@ namespace storm { template void ParameterRegionParser::parseParameterBoundaries(Valuation& lowerBoundaries, Valuation& upperBoundaries, std::string const& parameterBoundariesString, std::set const& consideredVariables) { - std::string::size_type positionOfFirstRelation = parameterBoundariesString.find("<="); STORM_LOG_THROW(positionOfFirstRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a '<=' after the first number"); std::string::size_type positionOfSecondRelation = parameterBoundariesString.find("<=", positionOfFirstRelation+2); STORM_LOG_THROW(positionOfSecondRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a '<=' after the parameter"); std::string parameter = parameterBoundariesString.substr(positionOfFirstRelation+2,positionOfSecondRelation-(positionOfFirstRelation+2)); + //removes all whitespaces from the parameter string: parameter.erase(std::remove_if (parameter.begin(), parameter.end(), ::isspace), parameter.end()); STORM_LOG_THROW(parameter.length()>0, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a parameter"); @@ -25,17 +26,19 @@ namespace storm { for (auto const& v : consideredVariables) { std::stringstream stream; stream << v; - std::string vAsString = stream.str(); if (parameter == stream.str()) { var = std::make_unique(v); + break; } } - STORM_LOG_ASSERT(var, "Could not find parameter " << parameter << " in the set of considered variables"); - - CoefficientType lb = storm::utility::convertNumber(parameterBoundariesString.substr(0,positionOfFirstRelation)); - CoefficientType ub = storm::utility::convertNumber(parameterBoundariesString.substr(positionOfSecondRelation+2)); - lowerBoundaries.emplace(std::make_pair(*var, lb)); - upperBoundaries.emplace(std::make_pair(*var, ub)); + if (var) { + CoefficientType lb = storm::utility::convertNumber(parameterBoundariesString.substr(0,positionOfFirstRelation)); + CoefficientType ub = storm::utility::convertNumber(parameterBoundariesString.substr(positionOfSecondRelation+2)); + lowerBoundaries.emplace(std::make_pair(*var, lb)); + upperBoundaries.emplace(std::make_pair(*var, ub)); + } else { + STORM_LOG_WARN("Could not find parameter " << parameter << " in the set of considered variables. Ignoring this parameter."); + } } template @@ -49,6 +52,12 @@ namespace storm { parseParameterBoundaries(lowerBoundaries, upperBoundaries, parameterBoundary, consideredVariables); } } + + // Check that all considered variables are bounded + for (auto const& v : consideredVariables) { + STORM_LOG_THROW(lowerBoundaries.count(v) > 0, storm::exceptions::WrongFormatException, "Variable " << v << " was not defined in region string."); + STORM_LOG_ASSERT(upperBoundaries.count(v) > 0, "Variable " << v << " has a lower but not an upper bound."); + } return storm::storage::ParameterRegion(std::move(lowerBoundaries), std::move(upperBoundaries)); } 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); diff --git a/src/storm/logic/Bound.h b/src/storm/logic/Bound.h index 21f2c7282..0b152ea06 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) const { + 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); }; 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..58970c77b 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; } @@ -205,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 91ccdf342..dc6396e71 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; @@ -109,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/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..322f36fdd --- /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 << (storm::solver::minimize(bv.first) ? "min" : "max") << " " << 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/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/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/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/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp index 07c732341..3a21d36db 100644 --- a/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp +++ b/src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp @@ -70,7 +70,7 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(env, pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); - STORM_LOG_THROW(pathFormula.getTimeBoundReference().isTimeBound(), storm::exceptions::NotImplementedException, "Currently step-bounded and reward=bpimded properties on MAs are not supported."); + STORM_LOG_THROW(pathFormula.getTimeBoundReference().isTimeBound(), storm::exceptions::NotImplementedException, "Currently step-bounded and reward-bounded properties on MAs are not supported."); double lowerBound = 0; double upperBound = 0; if (pathFormula.hasLowerBound()) { diff --git a/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp index f5137b0a7..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); @@ -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/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index b78ad3856..759d7a976 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.computeQuantile(env); + + 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/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 0e5cf1849..ff7cb5c98 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,17 +124,12 @@ 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; 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 0a0ea992c..3aff389fa 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,17 +122,12 @@ 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; 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 80059ce1e..6049785ec 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -3,21 +3,46 @@ #include #include "storm/storage/BitVector.h" +#include "storm/solver/OptimizationDirection.h" namespace storm { namespace modelchecker { 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 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; - 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; + + /// 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/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/EpochModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.cpp new file mode 100644 index 000000000..e8b059f32 --- /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 struct EpochModel; + template struct EpochModel; + template struct EpochModel; + template struct 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 efc3c81f7..3b1f36791 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -7,11 +7,14 @@ #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" #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 +34,9 @@ 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::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,23 +57,29 @@ namespace storm { objective.considersComplementaryEvent = false; objectives.push_back(std::move(objective)); - initialize(); + initialize(infinityBoundVariables); } template - void MultiDimensionalRewardUnfolding::initialize() { + 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); + + // collect which epoch steps are possible + possibleEpochSteps.clear(); + for (auto const& step : epochSteps) { + possibleEpochSteps.insert(step); + } } 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) { @@ -86,24 +96,43 @@ namespace storm { memLabel = "_" + memLabel; } dimension.memoryLabel = memLabel; - dimension.isUpperBounded = subformula.hasUpperBound(dim); - // for simplicity we do not allow intervals or 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()) { - dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); - dimension.scalingFactor = storm::utility::one(); + 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 + 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.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); + 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)); } @@ -112,8 +141,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.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(); @@ -150,7 +180,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); @@ -178,15 +208,9 @@ 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(); - + translateLowerBoundInfinityDimensions(epochSteps); } template @@ -201,7 +225,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(); @@ -216,29 +240,92 @@ 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()) { + // 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].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].boundType == DimensionBoundType::UpperBound)) { + discretizedBound -= storm::utility::one(); + } + } 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; + } + } + } + } + + 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().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. + 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; + } + } } - } 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; + 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() { + typename MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::getStartEpoch(bool setUnknownDimsToBottom) { 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].maxValue) { + epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } else { + 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); + } } STORM_LOG_TRACE("Start epoch is " << epochManager.toString(startEpoch)); return startEpoch; @@ -284,20 +371,26 @@ 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)); + 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())) { 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; } bool containsLowerBoundedObjective = false; for (auto const& dimension : dimensions) { - if (!dimension.isUpperBounded) { + if (dimension.boundType == DimensionBoundType::LowerBound) { containsLowerBoundedObjective = true; break; } @@ -327,7 +420,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; } @@ -409,7 +502,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); } } @@ -443,8 +536,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); @@ -544,9 +637,8 @@ 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; } - template template::type> @@ -605,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; @@ -728,7 +819,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) { @@ -755,11 +845,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); @@ -771,20 +859,19 @@ 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."); + 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]]; } 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())); } 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())); } @@ -798,6 +885,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 e1b96c64c..510cdcfd5 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 @@ -45,14 +36,27 @@ namespace storm { * */ MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::vector> const& objectives); - MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula); + + /*! + * 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; - 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); + EpochModel& setCurrentEpoch(Epoch const& epoch); void setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat); @@ -73,14 +77,21 @@ 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: void setCurrentEpochClass(Epoch const& epoch); - void initialize(); + 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); @@ -117,13 +128,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; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index af41c8a4f..aa387d084 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) { @@ -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); @@ -112,7 +111,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 +167,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 +234,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 +326,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 +392,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; } @@ -458,10 +481,15 @@ 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) { - 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].maxValue) { + epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get()); + } else { + epochManager.setBottomDimension(startEpoch, dim); + } } std::set seenEpochs({startEpoch}); @@ -481,28 +509,45 @@ namespace storm { } } } + + // 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].boundType != DimensionBoundType::Unbounded && !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 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].boundType != DimensionBoundType::Unbounded && dimensions[dim].maxValue) { + 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) { @@ -580,7 +625,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"); @@ -603,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 3f5b2a9a7..6ae408930 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h @@ -46,11 +46,13 @@ 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: - 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; @@ -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 new file mode 100644 index 000000000..3a5ca0239 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp @@ -0,0 +1,605 @@ +#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h" + +#include +#include +#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" +#include "storm/storage/expressions/Expressions.h" +#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" + +#include "storm/exceptions/NotSupportedException.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 + 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 + // Bounds are either constants or variables that are declared in the quantile formula. + // Prop op has optimality type + // No Prmin with lower cost bounds: Ec construction fails. In single obj we would do 1-Prmax[F "nonRewOrNonGoalEC"] but this invalidates other lower/upper cost bounds. + /* Todo: Current assumptions: + * Subformula is always prob operator with bounded until + * Each bound variable occurs at most once (change this?) + * cost bounds are assumed to be always non-strict (the epochs returned by the reward unfolding assume strict lower and non-strict upper cost bounds, though...) + * 'reasonable' quantile (e.g. not quantile(max B, Pmax>0.5 [F <=B G]) (we could also filter out these things early on) + */ + } + + 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(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; + + for (uint64_t dim = 0; dim < origBoundedUntil.getDimension(); ++dim) { + 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(dim)) { + lowerBounds.push_back(storm::logic::TimeBound(origBoundedUntil.isLowerBoundStrict(dim), origBoundedUntil.getLowerBound(dim))); + } else { + lowerBounds.push_back(boost::none); + } + if (origBoundedUntil.hasUpperBound(dim)) { + 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); + } 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); + } + if (transformations[dim] == BoundTransformation::LessEqualZero) { + lowerBounds.push_back(boost::none); + 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); + } + } + } + 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()); + } + + /// 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(); + } + + template + storm::storage::BitVector QuantileHelper::getOpenDimensions() const { + auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); + storm::storage::BitVector res(getDimension(), false); + for (uint64_t dim = 0; dim < getDimension(); ++dim) { + 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 { + storm::expressions::Variable const& dimVar = getVariableForDimension(dim); + 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::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 { + auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula(); + storm::storage::BitVector result(boundedUntil.getDimension(), false); + for (uint64_t dim = 0; dim < boundedUntil.getDimension(); ++dim) { + if (boundedUntil.hasLowerBound(dim) && boundedUntil.getLowerBound(dim).isVariable() && boundedUntil.getLowerBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) { + result.set(dim, true); + } + if (boundedUntil.hasUpperBound(dim) && boundedUntil.getUpperBound(dim).isVariable() && boundedUntil.getUpperBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) { + result.set(dim, true); + } + } + return result; + } + + template + std::vector> QuantileHelper::computeQuantile(Environment const& env) const { + std::vector> result; + Environment envCpy = env; // It might be necessary to increase the precision during the computation + numCheckedEpochs = 0; numPrecisionRefinements = 0; + if (getOpenDimensions().getNumberOfSetBits() == 1) { + uint64_t dimension = *getOpenDimensions().begin(); + + bool zeroSatisfiesFormula = checkLimitValue(envCpy, storm::storage::BitVector(getDimension(), false)); + bool infSatisfiesFormula = checkLimitValue(envCpy, getOpenDimensions()); + if (zeroSatisfiesFormula != infSatisfiesFormula) { + while (true) { + auto quantileRes = computeQuantileForDimension(envCpy, dimension); + if (quantileRes) { + result = {{storm::utility::convertNumber(quantileRes->first) * quantileRes->second}}; + break; + } + increasePrecision(envCpy); + ++numPrecisionRefinements; + } + } else if (zeroSatisfiesFormula) { + // thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula + if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) { + result = {{ storm::utility::zero() }}; + } else { + result = {{ storm::utility::infinity()}}; + } + } else { + // i.e., no bound value satisfies the formula + result = {{}}; + } + } else if (getOpenDimensions().getNumberOfSetBits() == 2) { + result = computeTwoDimensionalQuantile(envCpy); + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The quantile formula considers " << getOpenDimensions().getNumberOfSetBits() << " open dimensions. Only one- or two-dimensional quantiles are supported."); + } + 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; + } + + template + void filterDominatedPoints(std::vector>& points, std::vector const& dirs) { + std::vector> result; + // Note: this is slow and not inplace but also most likely not performance critical + for (auto const& p1 : points) { + bool p1Dominated = false; + for (auto const& p2 : points) { + assert(p1.size() == p2.size()); + bool p2DominatesP1 = false; + for (uint64_t i = 0; i < dirs.size(); ++i) { + if (storm::solver::minimize(dirs[i]) ? p2[i] <= p1[i] : p2[i] >= p1[i]) { + if (p2[i] != p1[i]) { + p2DominatesP1 = true; + } + } else { + p2DominatesP1 = false; + break; + } + } + if (p2DominatesP1) { + p1Dominated = true; + break; + } + } + if (!p1Dominated) { + result.push_back(p1); + } + } + points = std::move(result); + } + + template + std::vector> QuantileHelper::computeTwoDimensionalQuantile(Environment& env) const { + std::vector> result; + + auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); + auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula(); + std::vector dimensionsAsBitVector; + std::vector dimensions; + std::vector optimizationDirections; + std::vector lowerCostBounds; + for (auto const& dirVar : quantileFormula.getBoundVariables()) { + dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second)); + STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'."); + dimensions.push_back(*dimensionsAsBitVector.back().begin()); + lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back())); + optimizationDirections.push_back(dirVar.first); + } + STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions."); + if (optimizationDirections[0] == optimizationDirections[1]) { + if (lowerCostBounds[0] == lowerCostBounds[1]) { + // TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0]) + + bool infInfProbSatisfiesFormula = checkLimitValue(env, storm::storage::BitVector(getDimension(), false)); + bool zeroZeroProbSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1]); + if (infInfProbSatisfiesFormula != zeroZeroProbSatisfiesFormula) { + std::vector> singleQuantileValues; + for (uint64_t i = 0; i < 2; ++i) { + // find out whether the bounds B_i = 0 and B_1-i = infinity satisfy the formula + bool zeroInfSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[1-i]); + if (zeroInfSatisfiesFormula == storm::solver::minimize(optimizationDirections[0])) { + // There is bound value b such that the point B_i=0 and B_1-i = b is part of the result + singleQuantileValues.emplace_back(0, storm::utility::zero()); + } else { + // Compute quantile where 1-i is set to infinity + singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]).get()); + if (!storm::solver::minimize(optimizationDirections[i])) { + // When maximizing bounds, the computed single dimensional quantile value is sat for all values of the other bound. + // Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i] + ++singleQuantileValues.back().first; + } + } + // Decrease the value for lower cost bounds to convert >= to > + if (lowerCostBounds[i]) { + --singleQuantileValues[i].first; + } + } + std::vector currentEpochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0). + while (!exploreTwoDimensionalQuantile(env, singleQuantileValues, currentEpochValues, result)) { + increasePrecision(env); + ++numPrecisionRefinements; + } + } else if (infInfProbSatisfiesFormula) { + // then also zeroZeroProb satisfies the formula, i.e., all bound values satisfy the formula + if (storm::solver::minimize(optimizationDirections[0])) { + result = {{storm::utility::zero(), storm::utility::zero()}}; + } else { + result = {{ storm::utility::infinity(), storm::utility::infinity()}}; + } + } else { + // i.e., no bound value satisfies the formula + result = {{}}; + } + } else { + // TODO: this is an "unreasonable" case + } + } else { + // TODO: find reasonable min/max cases + } + return result; + } + + template + bool QuantileHelper::exploreTwoDimensionalQuantile(Environment const& env, std::vector> const& startEpochValues, std::vector& currentEpochValues, std::vector>& resultPoints) const { + // Initialize some data for easy access + auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula(); + auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula(); + std::vector dimensionsAsBitVector; + std::vector dimensions; + std::vector optimizationDirections; + std::vector lowerCostBounds; + for (auto const& dirVar : quantileFormula.getBoundVariables()) { + dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second)); + STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'."); + dimensions.push_back(*dimensionsAsBitVector.back().begin()); + lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back())); + optimizationDirections.push_back(dirVar.first); + } + STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions."); + + + MultiDimensionalRewardUnfolding rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector(getDimension(), BoundTransformation::None))); + + // initialize data that will be needed for each epoch + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + std::vector x, b; + std::unique_ptr> minMaxSolver; + + if (currentEpochValues[0] < 0 && currentEpochValues[1] < 0) { + // This case can only happen in these cases: + assert(lowerCostBounds[0]); + assert(currentEpochValues[0] == -1); + assert(currentEpochValues[1] == -1); + // This case has been checked already, so we can skip it. + // Skipping this is actually necessary, since the rewardUnfolding will handle formulas like [F{"a"}>A,{"b"}>B "init"] incorrectly if A=B=-1. + ++currentEpochValues[1]; + } + std::set exploredEpochs; + while (true) { + auto currEpoch = rewardUnfolding.getStartEpoch(true); + for (uint64_t i = 0; i < 2; ++i) { + if (currentEpochValues[i] >= 0) { + rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) currentEpochValues[i]); + } else { + rewardUnfolding.getEpochManager().setBottomDimension(currEpoch, dimensions[i]); + } + } + auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch); + for (auto const& epoch : epochOrder) { + if (exploredEpochs.count(epoch) == 0) { + auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); + ++numCheckedEpochs; + exploredEpochs.insert(epoch); + } + } + ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch); + bool propertySatisfied; + if (env.solver().isForceSoundness()) { + auto lowerUpperValue = getLowerUpperBound(env, currValue); + propertySatisfied = probOpFormula.getBound().isSatisfied(lowerUpperValue.first); + if (propertySatisfied != probOpFormula.getBound().isSatisfied(lowerUpperValue.second)) { + // unclear result due to insufficient precision. + return false; + } + } else { + propertySatisfied = probOpFormula.getBound().isSatisfied(currValue); + } + + // If the reward bounds should be as small as possible, we should stop as soon as the property is satisfied. + // If the reward bounds should be as large as possible, we should stop as soon as the property is violated and then go a step backwards + if (storm::solver::minimize(optimizationDirections[0]) == propertySatisfied) { + // We found another point for the result! Translate it to the original domain + auto point = currentEpochValues; + std::vector convertedPoint; + for (uint64_t i = 0; i < 2; ++i) { + if (lowerCostBounds[i]) { + // Translate > to >= + ++point[i]; + } + if (i == 1 && storm::solver::maximize(optimizationDirections[i])) { + // When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction + --point[i]; + } + if (point[i] < 0) { + // Skip this point + convertedPoint.clear(); + continue; + } + convertedPoint.push_back(storm::utility::convertNumber(point[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor); + } + if (!convertedPoint.empty()) { + resultPoints.push_back(std::move(convertedPoint)); + } + + if (currentEpochValues[1] == startEpochValues[1].first) { + break; + } else { + ++currentEpochValues[0]; + currentEpochValues[1] = startEpochValues[1].first; + } + } else { + ++currentEpochValues[1]; + } + } + + + // When maximizing, there are border cases where one dimension can be arbitrarily large + for (uint64_t i = 0; i < 2; ++i) { + if (storm::solver::maximize(optimizationDirections[i]) && startEpochValues[i].first > 0) { + std::vector newResultPoint(2); + newResultPoint[i] = storm::utility::convertNumber(startEpochValues[i].first - 1) * startEpochValues[i].second; + newResultPoint[1-i] = storm::utility::infinity(); + resultPoints.push_back(newResultPoint); + } + } + + filterDominatedPoints(resultPoints, optimizationDirections); + return true; + } + + template + boost::optional> QuantileHelper::computeQuantileForDimension(Environment const& env, uint64_t dimension) const { + // We assume that there is one bound value that violates the quantile and one bound value that satisfies it. + + // Let all other open bounds approach infinity + std::set infinityVariables; + for (auto const& d : getOpenDimensions()) { + if (d != dimension) { + infinityVariables.insert(getVariableForDimension(d)); + } + } + auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector(getDimension(), BoundTransformation::None)); + MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, infinityVariables); + + bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(dimension); + bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension)); + + // initialize data that will be needed for each epoch + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + std::vector x, b; + std::unique_ptr> minMaxSolver; + + uint64_t epochValue = 0; + std::set exploredEpochs; + while (true) { + auto currEpoch = rewardUnfolding.getStartEpoch(true); + rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimension, (uint64_t) epochValue); + auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch); + for (auto const& epoch : epochOrder) { + if (exploredEpochs.count(epoch) == 0) { + auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); + ++numCheckedEpochs; + exploredEpochs.insert(epoch); + } + } + ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch); + + bool propertySatisfied; + if (env.solver().isForceSoundness()) { + auto lowerUpperValue = getLowerUpperBound(env, currValue); + propertySatisfied = transformedFormula->getBound().isSatisfied(lowerUpperValue.first); + if (propertySatisfied != transformedFormula->getBound().isSatisfied(lowerUpperValue.second)) { + // unclear result due to insufficient precision. + return boost::none; + } + } else { + propertySatisfied = transformedFormula->getBound().isSatisfied(currValue); + } + // If the reward bound should be as small as possible, we should stop as soon as the property is satisfied. + // If the reward bound should be as large as possible, we should stop as soon as sthe property is violated and then go a step backwards + if (minimizingRewardBound && propertySatisfied) { + // We found a satisfying value! + if (!upperCostBound) { + // The rewardunfolding assumes strict lower cost bounds while we assume non-strict ones. Hence, >B becomes >=(B+1). + ++epochValue; + } + break; + } else if (!minimizingRewardBound && !propertySatisfied) { + // We found a non-satisfying value. Go one step back to get the largest satisfying value. + // ... however, lower cost bounds need to be converted from strict bounds >B to non strict bounds >=(B+1). + // Hence, no operation is necessary in this case. + if (upperCostBound) { + STORM_LOG_ASSERT(epochValue > 0, "The property does not seem to be satisfiable. This case should have been excluded earlier."); + --epochValue; + } + break; + } + ++epochValue; + } + return std::make_pair(epochValue, rewardUnfolding.getDimension(dimension).scalingFactor); + } + + template + bool QuantileHelper::checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const { + auto const& probabilityBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound(); + // Increase the precision until we get a conclusive result + while (true) { + ValueType numericResult = computeLimitValue(env, infDimensions); + if (env.solver().isForceSoundness()) { + auto lowerUpper = getLowerUpperBound(env, numericResult); + bool lowerSat = probabilityBound.isSatisfied(lowerUpper.first); + bool upperSat = probabilityBound.isSatisfied(lowerUpper.second); + if (lowerSat == upperSat) { + return lowerSat; + } else { + increasePrecision(env); + ++numPrecisionRefinements; + } + } else { + return probabilityBound.isSatisfied(numericResult); + } + } + } + + template + typename ModelType::ValueType QuantileHelper::computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const { + // To compute the limit for an upper bounded dimension, we can simply drop the bound + // To compute the limit for a lower bounded dimension, we only consider epoch steps that lie in an end component and check for the bound >0 instead. + // Notice, however, that this approach becomes problematic if, at the same time, we consider an upper bounded dimension with bound value zero. + std::vector bts(getDimension(), BoundTransformation::None); + std::set infinityVariables; + for (auto const& d : getOpenDimensions()) { + bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d); + if (infDimensions.get(d)) { + infinityVariables.insert(getVariableForDimension(d)); + } else { + if (upperCostBound) { + bts[d] = BoundTransformation::LessEqualZero; + } else { + bts[d] = BoundTransformation::GreaterEqualZero; + } + } + } + auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts); + + MultiDimensionalRewardUnfolding rewardUnfolding(model, transformedFormula, infinityVariables); + if (!rewardUnfolding.getProb1Objectives().empty()) { + assert(rewardUnfolding.getProb1Objectives().size() == 1); + // The probability is one. + return storm::utility::one(); + } + // Get lower and upper bounds for the solution. + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + + // Initialize epoch models + auto initEpoch = rewardUnfolding.getStartEpoch(); + auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch); + // initialize data that will be needed for each epoch + std::vector x, b; + std::unique_ptr> minMaxSolver; + + for (auto const& epoch : epochOrder) { + auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); + rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound)); + ++numCheckedEpochs; + } + + ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch); + STORM_LOG_TRACE("Limit probability for infinity dimensions " << infDimensions << " is " << numericResult << "."); + return numericResult; + } + + template class QuantileHelper>; + template class QuantileHelper>; + + } + } + } +} diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h new file mode 100644 index 000000000..13676a644 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h @@ -0,0 +1,70 @@ +#pragma once + +#include "storm/logic/QuantileFormula.h" +#include "boost/optional.hpp" + +namespace storm { + class Environment; + + namespace storage { + class BitVector; + } + + namespace modelchecker { + namespace helper { + namespace rewardbounded { + + template + class QuantileHelper { + typedef typename ModelType::ValueType ValueType; + public: + QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula); + + std::vector> computeQuantile(Environment const& env) const; + + private: + std::vector> computeTwoDimensionalQuantile(Environment& env) const; + 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. + /// 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 + */ + uint64_t getDimension() const; + + /*! + * 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; + storm::expressions::Variable const& getVariableForDimension(uint64_t const& dim) const; + + ModelType const& model; + storm::logic::QuantileFormula const& quantileFormula; + + /// Statistics + mutable uint64_t numCheckedEpochs; + mutable uint64_t numPrecisionRefinements; + }; + } + } + } +} 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; diff --git a/src/storm/models/sparse/Model.cpp b/src/storm/models/sparse/Model.cpp index 42b20db74..66da51770 100644 --- a/src/storm/models/sparse/Model.cpp +++ b/src/storm/models/sparse/Model.cpp @@ -40,9 +40,9 @@ namespace storm { // general components for all model types. STORM_LOG_THROW(this->getTransitionMatrix().getColumnCount() == stateCount, storm::exceptions::IllegalArgumentException, "Invalid column count of transition matrix."); STORM_LOG_ASSERT(components.rateTransitions || this->hasParameters() || this->getTransitionMatrix().isProbabilistic(), "The matrix is not probabilistic."); - STORM_LOG_THROW(this->getStateLabeling().getNumberOfItems() == stateCount, storm::exceptions::IllegalArgumentException, "Invalid item count of state labeling."); + STORM_LOG_THROW(this->getStateLabeling().getNumberOfItems() == stateCount, storm::exceptions::IllegalArgumentException, "Invalid item count (" << this->getStateLabeling().getNumberOfItems() << ") of state labeling (states: " << stateCount << ")."); for (auto const& rewardModel : this->getRewardModels()) { - STORM_LOG_THROW(!rewardModel.second.hasStateRewards() || rewardModel.second.getStateRewardVector().size() == stateCount, storm::exceptions::IllegalArgumentException, "Invalid size of state reward vector."); + STORM_LOG_THROW(!rewardModel.second.hasStateRewards() || rewardModel.second.getStateRewardVector().size() == stateCount, storm::exceptions::IllegalArgumentException, "Invalid size (" << rewardModel.second.getStateRewardVector().size() << ") of state reward vector (states:" << stateCount << ")."); STORM_LOG_THROW(!rewardModel.second.hasStateActionRewards() || rewardModel.second.getStateActionRewardVector().size() == choiceCount, storm::exceptions::IllegalArgumentException, "Invalid size of state reward vector."); STORM_LOG_ASSERT(!rewardModel.second.hasTransitionRewards() || rewardModel.second.getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix()), "The transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist."); } @@ -473,8 +473,15 @@ namespace storm { } return result; } + + std::set getAllParameters(Model const& model) { + std::set parameters = getProbabilityParameters(model); + std::set rewardParameters = getRewardParameters(model); + parameters.insert(rewardParameters.begin(), rewardParameters.end()); + return parameters; + } #endif - + template class Model; template class Model; diff --git a/src/storm/models/sparse/Model.h b/src/storm/models/sparse/Model.h index cc2fd3779..4e1d91b0a 100644 --- a/src/storm/models/sparse/Model.h +++ b/src/storm/models/sparse/Model.h @@ -419,10 +419,28 @@ namespace storm { boost::optional> choiceOrigins; }; - + #ifdef STORM_HAVE_CARL + /*! + * Get all probability parameters occurring on transitions. + * @param model Model. + * @return Set of parameters. + */ std::set getProbabilityParameters(Model const& model); + + /*! + * Get all parameters occurring in rewards. + * @param model Model. + * @return Set of parameters. + */ std::set getRewardParameters(Model const& model); + + /*! + * Get all parameters (probability and rewards) occurring in the model. + * @param model Model. + * @return Set of parameters. + */ + std::set getAllParameters(Model const& model); #endif } // namespace sparse } // namespace models diff --git a/src/storm/storage/MaximalEndComponentDecomposition.cpp b/src/storm/storage/MaximalEndComponentDecomposition.cpp index 7cbf36fd7..ec5484912 100644 --- a/src/storm/storage/MaximalEndComponentDecomposition.cpp +++ b/src/storm/storage/MaximalEndComponentDecomposition.cpp @@ -74,10 +74,10 @@ namespace storm { if (states) { endComponentStateSets.emplace_back(states->begin(), states->end(), true); } else { - std::vector states; - states.resize(transitionMatrix.getRowGroupCount()); - std::iota(states.begin(), states.end(), 0); - endComponentStateSets.emplace_back(states.begin(), states.end(), true); + std::vector allStates; + allStates.resize(transitionMatrix.getRowGroupCount()); + std::iota(allStates.begin(), allStates.end(), 0); + endComponentStateSets.emplace_back(allStates.begin(), allStates.end(), true); } storm::storage::BitVector statesToCheck(numberOfStates); storm::storage::BitVector includedChoices; diff --git a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp index 019ec40c2..2abe49b1e 100644 --- a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp +++ b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp @@ -211,7 +211,7 @@ namespace storm { } // Finally construct the quotient model. - this->quotient = std::make_shared(builder.build(), std::move(newLabeling), std::move(rewardModels)); + this->quotient = std::make_shared(builder.build(0,this->size(), this->size()), std::move(newLabeling), std::move(rewardModels)); } template 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; diff --git a/src/storm/transformer/EndComponentEliminator.h b/src/storm/transformer/EndComponentEliminator.h index 378edc86a..c522e1987 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); @@ -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; diff --git a/src/storm/utility/graph.cpp b/src/storm/utility/graph.cpp index d17d95edc..1515a31f7 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(state)) { + 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