Merge branch 'master' into dft

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

src/storm-pars/parser/ParameterRegionParser.cpp

@@ -1,3 +1,4 @@
+#include <storm/exceptions/WrongFormatException.h>
 #include "storm-pars/parser/ParameterRegionParser.h"
 
 #include "storm/utility/macros.h"
@@ -10,13 +11,13 @@ namespace storm {
 
         template<typename ParametricType>
         void ParameterRegionParser<ParametricType>::parseParameterBoundaries(Valuation& lowerBoundaries, Valuation& upperBoundaries, std::string const& parameterBoundariesString, std::set<VariableType> 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<VariableType>(v);
+                    break;
                 }
             }
-            STORM_LOG_ASSERT(var, "Could not find parameter " << parameter << " in the set of considered variables");
-            
-            CoefficientType lb = storm::utility::convertNumber<CoefficientType>(parameterBoundariesString.substr(0,positionOfFirstRelation));
-            CoefficientType ub = storm::utility::convertNumber<CoefficientType>(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<CoefficientType>(parameterBoundariesString.substr(0,positionOfFirstRelation));
+                CoefficientType ub = storm::utility::convertNumber<CoefficientType>(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<typename ParametricType>
@@ -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<ParametricType>(std::move(lowerBoundaries), std::move(upperBoundaries));
         }
 

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<storm::solver::OptimizationDirection, storm::expressions::Variable> 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<storm::solver::OptimizationDirection, storm::expressions::Variable> result(dir, var);
+            return result;
+        }
+
+        std::shared_ptr<storm::logic::Formula const> FormulaParserGrammar::createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula) {
+            return std::shared_ptr<storm::logic::Formula const>(new storm::logic::QuantileFormula(boundVariables, subformula));
+        }
+        
         std::set<storm::expressions::Variable> FormulaParserGrammar::getUndefinedConstants(std::shared_ptr<storm::logic::Formula const> const& formula) const {
             std::set<storm::expressions::Variable> result;
             std::set<storm::expressions::Variable> usedVariables = formula->getUsedVariables();

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<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> longRunAverageRewardFormula;
             
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> multiFormula;
+            qi::rule<Iterator, std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>(), qi::locals<storm::solver::OptimizationDirection>, Skipper> quantileBoundVariable;
+            qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> quantileFormula;
             
             // Parser that is used to recognize doubles only (as opposed to Spirit's double_ parser).
             boost::spirit::qi::real_parser<double, boost::spirit::qi::strict_real_policies<double>> strict_double;
@@ -229,6 +232,8 @@ namespace storm {
             std::shared_ptr<storm::logic::Formula const> createBinaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& leftSubformula, std::shared_ptr<storm::logic::Formula const> const& rightSubformula, storm::logic::BinaryBooleanStateFormula::OperatorType operatorType);
             std::shared_ptr<storm::logic::Formula const> createUnaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& subformula, boost::optional<storm::logic::UnaryBooleanStateFormula::OperatorType> const& operatorType);
             std::shared_ptr<storm::logic::Formula const> createMultiFormula(std::vector<std::shared_ptr<storm::logic::Formula const>> const& subformulas);
+            std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName);
+            std::shared_ptr<storm::logic::Formula const> createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula);
             
             std::set<storm::expressions::Variable> getUndefinedConstants(std::shared_ptr<storm::logic::Formula const> const& formula) const;
             storm::jani::Property createProperty(boost::optional<std::string> const& propertyName, storm::modelchecker::FilterType const& filterType, std::shared_ptr<storm::logic::Formula const> const& formula, std::shared_ptr<storm::logic::Formula const> const& states);

src/storm/logic/Bound.h

@@ -16,6 +16,21 @@ namespace storm {
             ComparisonType comparisonType;
             storm::expressions::Expression threshold;
 
+            template <typename ValueType>
+            bool isSatisfied(ValueType const& compareValue) const {
+                ValueType thresholdAsValueType = storm::utility::convertNumber<ValueType>(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);
         };
 

src/storm/logic/CloneVisitor.cpp

@@ -108,6 +108,11 @@ namespace storm {
             return std::static_pointer_cast<Formula>(std::make_shared<MultiObjectiveFormula>(subformulas));
         }
         
+        boost::any CloneVisitor::visit(QuantileFormula const& f, boost::any const& data) const {
+            std::shared_ptr<Formula> subformula = boost::any_cast<std::shared_ptr<Formula>>(f.getSubformula().accept(*this, data));
+            return std::static_pointer_cast<Formula>(std::make_shared<QuantileFormula>(f.getBoundVariables(), subformula));
+        }
+        
         boost::any CloneVisitor::visit(NextFormula const& f, boost::any const& data) const {
             std::shared_ptr<Formula> subformula = boost::any_cast<std::shared_ptr<Formula>>(f.getSubformula().accept(*this, data));
             return std::static_pointer_cast<Formula>(std::make_shared<NextFormula>(subformula));

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;

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<MultiObjectiveFormula const&>(*this);
         }
         
+        QuantileFormula& Formula::asQuantileFormula() {
+            return dynamic_cast<QuantileFormula&>(*this);
+        }
+        
+        QuantileFormula const& Formula::asQuantileFormula() const {
+            return dynamic_cast<QuantileFormula const&>(*this);
+        }
+        
         BinaryStateFormula& Formula::asBinaryStateFormula() {
             return dynamic_cast<BinaryStateFormula&>(*this);
         }

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;
             

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<FormulaInformation>(f.getSubformula().accept(*this, data)).setContainsNextFormula();
         }

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;

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;

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"

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;

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<InheritedInformation const&>(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<InheritedInformation const&>(data);
             bool result = inherited.getSpecification().areNextFormulasAllowed();

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;

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;
         }

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();
 
     }
 }

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<bool>(f.getSubformula().accept(*this, data));
         }

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;

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<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> 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<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> 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<std::shared_ptr<AtomicExpressionFormula const>>& atomicExpressionFormulas) const {
+            subformula->gatherAtomicExpressionFormulas(atomicExpressionFormulas);
+        }
+        
+        void QuantileFormula::gatherAtomicLabelFormulas(std::vector<std::shared_ptr<AtomicLabelFormula const>>& atomicLabelFormulas) const {
+            subformula->gatherAtomicLabelFormulas(atomicLabelFormulas);
+        }
+        
+        void QuantileFormula::gatherReferencedRewardModels(std::set<std::string>& 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;
+        }
+    }
+}

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<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> 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<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> 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<std::shared_ptr<AtomicExpressionFormula const>>& atomicExpressionFormulas) const override;
+            virtual void gatherAtomicLabelFormulas(std::vector<std::shared_ptr<AtomicLabelFormula const>>& atomicLabelFormulas) const override;
+            virtual void gatherReferencedRewardModels(std::set<std::string>& referencedRewardModels) const override;
+            
+            virtual std::ostream& writeToStream(std::ostream& out) const override;
+        private:
+            std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> boundVariables;
+            std::shared_ptr<Formula const> subformula;
+        };
+    }
+}

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;
             

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;
             }

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.");
         }

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;

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<typename ModelType>
+        std::unique_ptr<CheckResult> AbstractModelChecker<ModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> 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.
         ///////////////////////////////////////////////

src/storm/modelchecker/AbstractModelChecker.h

@@ -91,7 +91,10 @@ namespace storm {
   
             // The methods to check multi-objective formulas.
             virtual std::unique_ptr<CheckResult> checkMultiObjectiveFormula(Environment const& env, CheckTask<storm::logic::MultiObjectiveFormula, ValueType> const& checkTask);
-                  
+  
+            // The methods to check quantile formulas.
+            virtual std::unique_ptr<CheckResult> checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask);
+            
         };
     }
 }

src/storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp

@@ -70,7 +70,7 @@ namespace storm {
             std::unique_ptr<CheckResult> 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()) {

src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.cpp

@@ -88,7 +88,7 @@ namespace storm {
                     if (storm::settings::getModule<storm::settings::modules::IOSettings>().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) {
                         std::vector<ValueType> 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<ValueType>(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor);
                         }
                         auto const& solution = rewardUnfolding.getInitialStateResult(epoch);
@@ -309,7 +309,7 @@ namespace storm {
             }
 
             template <class SparseMdpModelType>
-            void RewardBoundedMdpPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector) {
+            void RewardBoundedMdpPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(Environment const& env, helper::rewardbounded::EpochModel<ValueType, false> const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector) {
                 if (epochModel.epochMatrixChanged) {
                 
                     // Update the cached MinMaxSolver data

src/storm/modelchecker/multiobjective/pcaa/RewardBoundedMdpPcaaWeightVectorChecker.h

@@ -64,7 +64,7 @@ namespace storm {
                 
                 void computeEpochSolution(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData);
                 
-                void updateCachedData(Environment const& env, typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector);
+                void updateCachedData(Environment const& env, typename helper::rewardbounded::EpochModel<ValueType, false> const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector);
                 
                 storm::utility::Stopwatch swAll, swEpochModelBuild, swEpochModelAnalysis;
                 uint64_t numCheckedEpochs, numChecks;

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<typename SparseMdpModelType>
@@ -220,6 +226,22 @@ namespace storm {
             return multiobjective::performMultiObjectiveModelChecking(env, this->getModel(), checkTask.getFormula());
         }
         
+        template<typename SparseMdpModelType>
+        std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> 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<SparseMdpModelType> qHelper(this->getModel(), checkTask.getFormula());
+            auto res = qHelper.computeQuantile(env);
+            
+            if (res.size() == 1 && res.front().size() == 1) {
+                return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(initialState, std::move(res.front().front())));
+            } else {
+                return std::unique_ptr<CheckResult>(new ExplicitParetoCurveCheckResult<ValueType>(initialState, std::move(res)));
+            }
+        }
+        
         template class SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<double>>;
 
 #ifdef STORM_HAVE_CARL

src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h

@@ -33,6 +33,7 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeLongRunAverageProbabilities(Environment const& env, CheckTask<storm::logic::StateFormula, ValueType> const& checkTask) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
             virtual std::unique_ptr<CheckResult> checkMultiObjectiveFormula(Environment const& env, CheckTask<storm::logic::MultiObjectiveFormula, ValueType> const& checkTask) override;
+            virtual std::unique_ptr<CheckResult> checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) override;
             
         };
     } // namespace modelchecker

src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp

@@ -81,75 +81,6 @@ namespace storm {
                 return result;
             }
             
-            template<typename ValueType>
-            std::vector<ValueType> analyzeTrivialDtmcEpochModel(typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel) {
-                
-                std::vector<ValueType> 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<ValueType>());
-                        }
-                    }
-                }
-                return epochResult;
-            }
-            
-            template<typename ValueType>
-            std::vector<ValueType> analyzeNonTrivialDtmcEpochModel(Environment const& env, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>& linEqSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
- 
-                // Update some data for the case that the Matrix has changed
-                if (epochModel.epochMatrixChanged) {
-                    x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
-                    storm::solver::GeneralLinearEquationSolverFactory<ValueType> 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<ValueType>());
-                std::vector<ValueType> 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<storm::storage::sparse::state_type, storm::RationalFunction> SparseDtmcPrctlHelper<storm::RationalFunction>::computeRewardBoundedValues(Environment const& env, storm::models::sparse::Dtmc<storm::RationalFunction> const& model, std::shared_ptr<storm::logic::OperatorFormula const> 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<ValueType> 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<ValueType>(epochModel));
-                    } else {
-                        rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialDtmcEpochModel<ValueType>(preciseEnv, epochModel, x, b, linEqSolver, lowerBound, upperBound));
-                    }
+                    rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(preciseEnv, x, b, linEqSolver, lowerBound, upperBound));
                     swCheck.stop();
                     if (storm::settings::getModule<storm::settings::modules::IOSettings>().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) {
                         std::vector<ValueType> 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<ValueType>(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor);
                         }
                         cdfEntry.push_back(rewardUnfolding.getInitialStateResult(epoch));

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -90,103 +90,7 @@ namespace storm {
                 return result;
             }
 
-            template<typename ValueType>
-            std::vector<ValueType> analyzeTrivialMdpEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel) {
-                // Assert that the epoch model is indeed trivial
-                assert(epochModel.epochMatrix.getEntryCount() == 0);
-                
-                std::vector<ValueType> 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<ValueType>();
-                        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<typename ValueType>
-            std::vector<ValueType> analyzeNonTrivialMdpEpochModel(Environment const& env, OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
- 
-                // Update some data for the case that the Matrix has changed
-                if (epochModel.epochMatrixChanged) {
-                    x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
-                    storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> 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<ValueType>());
-                std::vector<ValueType> 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<typename ValueType>
             std::map<storm::storage::sparse::state_type, ValueType> SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(Environment const& env, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& 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<ValueType>(dir, epochModel));
-                    } else {
-                        rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialMdpEpochModel<ValueType>(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<storm::settings::modules::IOSettings>().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) {
                         std::vector<ValueType> 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<ValueType>(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor);
                         }
                         cdfEntry.push_back(rewardUnfolding.getInitialStateResult(epoch));

src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h

@@ -3,21 +3,46 @@
 #include <boost/optional.hpp>
 
 #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<typename ValueType>
                 struct Dimension {
+                    /// The formula describing this dimension
                     std::shared_ptr<storm::logic::Formula const> 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<std::string> 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<uint64_t> maxValue;
+
+                    /// Whether we minimize/maximize the objective for this dimension
+                    boost::optional<storm::solver::OptimizationDirection> optimizationDirection;
                 };
             }
         }

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;

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;

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<typename ValueType>
+                std::vector<ValueType> analyzeTrivialDtmcEpochModel(EpochModel<ValueType, true>& epochModel) {
+
+                    std::vector<ValueType> 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<ValueType>());
+                            }
+                        }
+                    }
+                    return epochResult;
+                }
+
+                template<typename ValueType>
+                std::vector<ValueType> analyzeNonTrivialDtmcEpochModel(Environment const& env, EpochModel<ValueType, true>& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>& linEqSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
+
+                    // Update some data for the case that the Matrix has changed
+                    if (epochModel.epochMatrixChanged) {
+                        x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
+                        storm::solver::GeneralLinearEquationSolverFactory<ValueType> 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<ValueType>());
+                    std::vector<ValueType> 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<typename ValueType>
+                std::vector<ValueType> analyzeTrivialMdpEpochModel(OptimizationDirection dir, EpochModel<ValueType, true>& epochModel) {
+                    // Assert that the epoch model is indeed trivial
+                    assert(epochModel.epochMatrix.getEntryCount() == 0);
+
+                    std::vector<ValueType> 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<ValueType>();
+                            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<typename ValueType>
+                std::vector<ValueType> analyzeNonTrivialMdpEpochModel(Environment const& env, OptimizationDirection dir, EpochModel<ValueType, true>& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
+
+                    // Update some data for the case that the Matrix has changed
+                    if (epochModel.epochMatrixChanged) {
+                        x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
+                        storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> 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<ValueType>());
+                    std::vector<ValueType> 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<double> EpochModel<double, true>::analyzeSingleObjective(
+                        const storm::Environment &env, std::vector<double> &x, std::vector<double> &b,
+                        std::unique_ptr<storm::solver::LinearEquationSolver<double>> &linEqSolver,
+                        const boost::optional<double> &lowerBound, const boost::optional<double> &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<double>(*this);
+                    } else {
+                        return analyzeNonTrivialDtmcEpochModel<double>(env, *this, x, b, linEqSolver, lowerBound, upperBound);
+                    }
+                }
+
+                template<>
+                std::vector<double> EpochModel<double, true>::analyzeSingleObjective(
+                        const storm::Environment &env, storm::OptimizationDirection dir, std::vector<double> &x,
+                        std::vector<double> &b,
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> &minMaxSolver,
+                        const boost::optional<double> &lowerBound, const boost::optional<double> &upperBound) {
+                    // If the epoch matrix is empty we do not need to solve a linear equation system
+                    if (epochMatrix.getEntryCount() == 0) {
+                        return analyzeTrivialMdpEpochModel<double>(dir, *this);
+                    } else {
+                        return analyzeNonTrivialMdpEpochModel<double>(env, dir, *this, x, b, minMaxSolver, lowerBound, upperBound);
+                    }
+                }
+
+                template<>
+                std::vector<storm::RationalNumber> EpochModel<storm::RationalNumber, true>::analyzeSingleObjective(
+                        const storm::Environment &env, std::vector<storm::RationalNumber> &x, std::vector<storm::RationalNumber> &b,
+                        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> &linEqSolver,
+                        const boost::optional<storm::RationalNumber> &lowerBound, const boost::optional<storm::RationalNumber> &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<storm::RationalNumber>(*this);
+                    } else {
+                        return analyzeNonTrivialDtmcEpochModel<storm::RationalNumber>(env, *this, x, b, linEqSolver, lowerBound, upperBound);
+                    }
+                }
+
+                template<>
+                std::vector<storm::RationalNumber> EpochModel<storm::RationalNumber, true>::analyzeSingleObjective(
+                        const storm::Environment &env, storm::OptimizationDirection dir, std::vector<storm::RationalNumber> &x,
+                        std::vector<storm::RationalNumber> &b,
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<storm::RationalNumber>> &minMaxSolver,
+                        const boost::optional<storm::RationalNumber> &lowerBound, const boost::optional<storm::RationalNumber> &upperBound) {
+                    // If the epoch matrix is empty we do not need to solve a linear equation system
+                    if (epochMatrix.getEntryCount() == 0) {
+                        return analyzeTrivialMdpEpochModel<storm::RationalNumber>(dir, *this);
+                    } else {
+                        return analyzeNonTrivialMdpEpochModel<storm::RationalNumber>(env, dir, *this, x, b, minMaxSolver, lowerBound, upperBound);
+                    }
+                }
+
+                template struct EpochModel<double, true>;
+                template struct EpochModel<double, false>;
+                template struct EpochModel<storm::RationalNumber, true>;
+                template struct EpochModel<storm::RationalNumber, false>;
+            }
+        }
+    }
+}

src/storm/modelchecker/prctl/helper/rewardbounded/EpochModel.h

@@ -0,0 +1,48 @@
+#pragma once
+
+#include <vector>
+#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<typename ValueType, bool SingleObjectiveMode>
+                struct EpochModel {
+                    typedef typename std::conditional<SingleObjectiveMode, ValueType, std::vector < ValueType>>::type SolutionType;
+
+                    bool epochMatrixChanged;
+                    storm::storage::SparseMatrix<ValueType> epochMatrix;
+                    storm::storage::BitVector stepChoices;
+                    std::vector<SolutionType> stepSolutions;
+                    std::vector<std::vector<ValueType>> objectiveRewards;
+                    std::vector<storm::storage::BitVector> 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<storm::solver::LinearEquationSolverProblemFormat> equationSolverProblemFormat;
+
+                    /*!
+                     * Analyzes the epoch model, i.e., solves the represented equation system. This method assumes a nondeterministic model.
+                     */
+                    std::vector<ValueType> analyzeSingleObjective(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound);
+
+                    /*!
+                     * Analyzes the epoch model, i.e., solves the represented equation system. This method assumes a deterministic model.
+                     */
+                    std::vector<ValueType> analyzeSingleObjective(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>& linEqSolver, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound);
+                };
+
+
+
+
+            }
+        }
+    }
+}

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<typename ValueType, bool SingleObjectiveMode>
-                MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::MultiDimensionalRewardUnfolding(storm::models::sparse::Model<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> objectiveFormula) : model(model) {
-                    
+                MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::MultiDimensionalRewardUnfolding(storm::models::sparse::Model<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> objectiveFormula, std::set<storm::expressions::Variable> 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<typename ValueType, bool SingleObjectiveMode>
-                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initialize() {
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initialize(std::set<storm::expressions::Variable> const& infinityBoundVariables) {
                     
                     maxSolutionsStored = 0;
                     
                     STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives.");
                     std::vector<Epoch> 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<typename ValueType, bool SingleObjectiveMode>
-                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initializeObjectives(std::vector<Epoch>& epochSteps) {
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initializeObjectives(std::vector<Epoch>& epochSteps, std::set<storm::expressions::Variable> const& infinityBoundVariables) {
                     std::vector<std::vector<uint64_t>> 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<uint64_t>(model.getTransitionMatrix().getRowCount(), 1));
-                                    dimension.scalingFactor = storm::utility::one<ValueType>();
+                                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<uint64_t>(model.getTransitionMatrix().getRowCount(), 0));
+                                    dimension.scalingFactor = storm::utility::zero<ValueType>();
+                                    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<ValueType> actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix());
-                                    auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector<ValueType, uint64_t>(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<uint64_t>(model.getTransitionMatrix().getRowCount(), 1));
+                                        dimension.scalingFactor = storm::utility::one<ValueType>();
+                                    } 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<ValueType> actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix());
+                                        auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector<ValueType, uint64_t>(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<ValueType> 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<uint64_t>(model.getTransitionMatrix().getRowCount(), 1));
                                     dimension.scalingFactor = storm::utility::one<ValueType>();
@@ -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<typename ValueType, bool SingleObjectiveMode>
@@ -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<ValueType>(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<ValueType>();
+                        
+                        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<ValueType>(bound.evaluateAsRational());
+                                discretizedBound /= dimensions[dim].scalingFactor;
+                                if (storm::utility::isInteger(discretizedBound)) {
+                                    if (isStrict == (dimensions[dim].boundType == DimensionBoundType::UpperBound)) {
+                                        discretizedBound -= storm::utility::one<ValueType>();
+                                    }
+                                } else {
+                                    discretizedBound = storm::utility::floor(discretizedBound);
+                                }
+                                uint64_t dimensionValue = storm::utility::convertNumber<uint64_t>(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<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::translateLowerBoundInfinityDimensions(std::vector<Epoch>& 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<uint64_t>(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<ValueType> 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 ValueType, bool SingleObjectiveMode>
-                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStartEpoch() {
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::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 ValueType, bool SingleObjectiveMode>
-                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::EpochModel& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpoch(Epoch const& epoch) {
+                EpochModel<ValueType, SingleObjectiveMode>& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpoch(Epoch const& epoch) {
                     STORM_LOG_DEBUG("Setting model for epoch " << epochManager.toString(epoch));
+                    STORM_LOG_THROW(getProb1Objectives().empty(), storm::exceptions::InvalidOperationException, "The objective " << objectives[*getProb1Objectives().begin()].formula << " is satisfied already in the initial state(s). This special case can not be handled by the reward unfolding."); // This case should have been handled 'from the outside'
                     
                     // Check if we need to update the current epoch class
                     if (!currentEpoch || !epochManager.compareEpochClass(epoch, currentEpoch.get())) {
                         setCurrentEpochClass(epoch);
                         epochModel.epochMatrixChanged = true;
+                        if (storm::settings::getModule<storm::settings::modules::CoreSettings>().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<typename ValueType, bool SingleObjectiveMode>
                 void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::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<typename ValueType, bool SingleObjectiveMode>
                 template<bool SO, typename std::enable_if<SO, int>::type>
@@ -605,7 +697,6 @@ namespace storm {
                 
                 template<typename ValueType, bool SingleObjectiveMode>
                 ValueType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) {
-                    // if (storm::settings::getModule<storm::settings::modules::GeneralSettings>().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<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::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 ValueType, bool SingleObjectiveMode>
                 typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::EpochSolution const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochSolution(std::map<Epoch, EpochSolution const*> const& solutions, Epoch const& epoch) {
                     auto epochSolutionIt = solutions.find(epoch);
@@ -771,20 +859,19 @@ namespace storm {
                 typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState) {
                     STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution at an unexisting product state.");
                     STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at a state for which no solution was stored.");
+                    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 ValueType, bool SingleObjectiveMode>
                 typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::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 ValueType, bool SingleObjectiveMode>
                 typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::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<typename ValueType, bool SingleObjectiveMode>
+                storm::storage::BitVector const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getProb1Objectives() const {
+                    return productModel->getProb1Objectives();
+                }
+                
                 template class MultiDimensionalRewardUnfolding<double, true>;
                 template class MultiDimensionalRewardUnfolding<double, false>;
                 template class MultiDimensionalRewardUnfolding<storm::RationalNumber, true>;

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<SingleObjectiveMode, ValueType, std::vector<ValueType>>::type SolutionType;
-    
-                    struct EpochModel {
-                        bool epochMatrixChanged;
-                        storm::storage::SparseMatrix<ValueType> epochMatrix;
-                        storm::storage::BitVector stepChoices;
-                        std::vector<SolutionType> stepSolutions;
-                        std::vector<std::vector<ValueType>> objectiveRewards;
-                        std::vector<storm::storage::BitVector> objectiveRewardFilter;
-                        storm::storage::BitVector epochInStates;
-                    };
-                    
+
                     /*
                      *
                      * @param model The (preprocessed) model
@@ -45,14 +36,27 @@ namespace storm {
                      *
                      */
                     MultiDimensionalRewardUnfolding(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
-                    MultiDimensionalRewardUnfolding(storm::models::sparse::Model<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> 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<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> objectiveFormula, std::set<storm::expressions::Variable> 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<Epoch> getEpochComputationOrder(Epoch const& startEpoch);
                     
-                    EpochModel& setCurrentEpoch(Epoch const& epoch);
+                    EpochModel<ValueType, SingleObjectiveMode>& setCurrentEpoch(Epoch const& epoch);
                     
                     void setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat);
                     
@@ -73,14 +77,21 @@ namespace storm {
                     
                     EpochManager const& getEpochManager() const;
                     Dimension<ValueType> const& getDimension(uint64_t dim) const;
+
+                    /*!
+                     * Returns objectives that are always satisfied (i.e., have probability one) in all initial states.
+                     * These objectives can not be handled by this as they can not be translated into expected rewards.
+                     */
+                    storm::storage::BitVector const& getProb1Objectives() const;
                     
                 private:
                 
                     void setCurrentEpochClass(Epoch const& epoch);
-                    void initialize();
+                    void initialize(std::set<storm::expressions::Variable> const& infinityBoundVariables = {});
                     
-                    void initializeObjectives(std::vector<Epoch>& epochSteps);
+                    void initializeObjectives(std::vector<Epoch>& epochSteps, std::set<storm::expressions::Variable> const& infinityBoundVariables);
                     void computeMaxDimensionValues();
+                    void translateLowerBoundInfinityDimensions(std::vector<Epoch>& epochSteps);
                     
                     void initializeMemoryProduct(std::vector<Epoch> const& epochSteps);
                     
@@ -117,13 +128,10 @@ namespace storm {
                     std::vector<uint64_t> epochModelToProductChoiceMap;
                     std::shared_ptr<std::vector<uint64_t> const> productStateToEpochModelInStateMap;
                     std::set<Epoch> possibleEpochSteps;
-                    
-                    EpochModel epochModel;
+
+                    EpochModel<ValueType, SingleObjectiveMode> epochModel;
                     boost::optional<Epoch> currentEpoch;
-                    
-                    // In case of DTMCs we have different options for the equation problem format the epoch model will have.
-                    boost::optional<storm::solver::LinearEquationSolverProblemFormat> equationSolverProblemFormatForEpochModel;
-                    
+
                     EpochManager epochManager;
                     
                     std::vector<Dimension<ValueType>> dimensions;

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

@@ -21,7 +21,7 @@ namespace storm {
             namespace rewardbounded {
                 
                 template<typename ValueType>
-                ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()) {
+                ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> 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<typename ValueType>
-                storm::storage::MemoryStructure ProductModel<ValueType>::computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const {
+                storm::storage::MemoryStructure ProductModel<ValueType>::computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) {
                     
                     storm::modelchecker::SparsePropositionalModelChecker<storm::models::sparse::Model<ValueType>> mc(model);
                     
@@ -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<EpochClass> initEpochClasses;
+                    initEpochClasses.push_back(epochManager.getEpochClass(epochManager.getZeroEpoch()));
+                    for (uint64_t dim = 0; dim < dimensions.size(); ++dim) {
+                        Dimension<ValueType> 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<EpochClass> 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<Epoch> possibleSteps(originalModelSteps.begin(), originalModelSteps.end());
@@ -303,7 +326,7 @@ namespace storm {
                 
                 template<typename ValueType>
                 uint64_t ProductModel<ValueType>::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<typename ValueType>
                 void ProductModel<ValueType>::collectReachableEpochClasses(std::set<EpochClass, std::function<bool(EpochClass const&, EpochClass const&)>>& reachableEpochClasses, std::set<Epoch> 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<Epoch> 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<EpochClass> 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<typename ValueType>
                 void ProductModel<ValueType>::computeReachableStates(EpochClass const& epochClass, std::vector<EpochClass> 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<typename ValueType>
+                storm::storage::BitVector const& ProductModel<ValueType>::getProb1Objectives() {
+                    return prob1Objectives;
+                }
+                
                 template class ProductModel<double>;
                 template class ProductModel<storm::RationalNumber>;
             }

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<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const;
+                    storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
                     std::vector<MemoryState> computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const;
     
                     
@@ -75,7 +77,7 @@ namespace storm {
                     std::vector<uint64_t> productToModelStateMap;
                     std::vector<MemoryState> productToMemoryStateMap;
                     std::vector<uint64_t> choiceToStateMap;
-                    
+                    storm::storage::BitVector prob1Objectives; /// Objectives that are already satisfied in the initial state
                 };
             }
         }

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

@@ -0,0 +1,605 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h"
+
+#include <set>
+#include <vector>
+#include <memory>
+#include <boost/optional.hpp>
+
+#include "storm/environment/solver/MinMaxSolverEnvironment.h"
+
+#include "storm/models/sparse/Mdp.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
+#include "storm/storage/expressions/ExpressionManager.h"
+#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<typename ModelType>
+                QuantileHelper<ModelType>::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<storm::logic::ProbabilityOperatorFormula> transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector<BoundTransformation> 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<std::shared_ptr<storm::logic::Formula const>> leftSubformulas, rightSubformulas;
+                    std::vector<boost::optional<storm::logic::TimeBound>> lowerBounds, upperBounds;
+                    std::vector<storm::logic::TimeBoundReference> 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<storm::logic::Formula> newBoundedUntil;
+                    if (origBoundedUntil.hasMultiDimensionalSubformulas()) {
+                        newBoundedUntil = std::make_shared<storm::logic::BoundedUntilFormula>(leftSubformulas, rightSubformulas, lowerBounds, upperBounds, timeBoundReferences);
+                    } else {
+                        newBoundedUntil = std::make_shared<storm::logic::BoundedUntilFormula>(origBoundedUntil.getLeftSubformula().asSharedPointer(), origBoundedUntil.getRightSubformula().asSharedPointer(), lowerBounds, upperBounds, timeBoundReferences);
+                    }
+                    return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, boundedUntilOperator.getOperatorInformation());
+                }
+
+                /// Increases the precision of solver results
+                void increasePrecision(storm::Environment& env) {
+                    STORM_LOG_DEBUG("Increasing precision of underlying solver.");
+                    auto factor = storm::utility::convertNumber<storm::RationalNumber, std::string>("0.1");
+                    env.solver().setLinearEquationSolverPrecision(static_cast<storm::RationalNumber>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get() * factor));
+                    env.solver().minMax().setPrecision(env.solver().minMax().getPrecision() * factor);
+                }
+                
+                /// Computes a lower/ upper bound on the actual result of a minmax or linear equation solver
+                template<typename ValueType>
+                std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& value, bool minMax = true) {
+                    ValueType prec;
+                    bool relative;
+                    if (minMax) {
+                        prec = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision());
+                        relative = env.solver().minMax().getRelativeTerminationCriterion();
+                    } else {
+                        prec = storm::utility::convertNumber<ValueType>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get());
+                        relative = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).second.get();
+                    }
+                    if (relative) {
+                        return std::make_pair<ValueType, ValueType>(value * (1/(prec + 1)), value * (1 + prec/(prec +1)));
+                    } else {
+                        return std::make_pair<ValueType, ValueType>(value - prec, value + prec);
+                    }
+                }
+                
+                template<typename ModelType>
+                uint64_t QuantileHelper<ModelType>::getDimension() const {
+                    return quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().getDimension();
+                }
+
+                template<typename ModelType>
+                storm::storage::BitVector QuantileHelper<ModelType>::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<typename ModelType>
+                storm::solver::OptimizationDirection const& QuantileHelper<ModelType>::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<typename ModelType>
+                storm::expressions::Variable const& QuantileHelper<ModelType>::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<typename ModelType>
+                storm::storage::BitVector QuantileHelper<ModelType>::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<typename ModelType>
+                std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) const {
+                    std::vector<std::vector<ValueType>> 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<ValueType>(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<ValueType>() }};
+                            } else {
+                                result = {{ storm::utility::infinity<ValueType>()}};
+                            }
+                        } 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<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
+                        std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl;
+                        std::cout << "Number of required precision refinements: " << numPrecisionRefinements << std::endl;
+                    }
+                    return result;
+                }
+
+                template<typename ValueType>
+                void filterDominatedPoints(std::vector<std::vector<ValueType>>& points, std::vector<storm::solver::OptimizationDirection> const& dirs) {
+                    std::vector<std::vector<ValueType>> 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<typename ModelType>
+                std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment& env) const {
+                    std::vector<std::vector<ValueType>> result;
+
+                    auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
+                    auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
+                    std::vector<storm::storage::BitVector> dimensionsAsBitVector;
+                    std::vector<uint64_t> dimensions;
+                    std::vector<storm::solver::OptimizationDirection> optimizationDirections;
+                    std::vector<bool> lowerCostBounds;
+                    for (auto const& dirVar : quantileFormula.getBoundVariables()) {
+                        dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
+                        STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
+                        dimensions.push_back(*dimensionsAsBitVector.back().begin());
+                        lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
+                        optimizationDirections.push_back(dirVar.first);
+                    }
+                    STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
+                    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<std::pair<int64_t, typename ModelType::ValueType>> 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<ValueType>());
+                                    } 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<int64_t> currentEpochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0).
+                                while (!exploreTwoDimensionalQuantile(env, singleQuantileValues, currentEpochValues, result)) {
+                                    increasePrecision(env);
+                                    ++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<ValueType>(), storm::utility::zero<ValueType>()}};
+                                } else {
+                                    result = {{ storm::utility::infinity<ValueType>(),  storm::utility::infinity<ValueType>()}};
+                                }
+                            } 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<typename ModelType>
+                bool QuantileHelper<ModelType>::exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const {
+                    // Initialize some data for easy access
+                    auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
+                    auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
+                    std::vector<storm::storage::BitVector> dimensionsAsBitVector;
+                    std::vector<uint64_t> dimensions;
+                    std::vector<storm::solver::OptimizationDirection> optimizationDirections;
+                    std::vector<bool> lowerCostBounds;
+                    for (auto const& dirVar : quantileFormula.getBoundVariables()) {
+                        dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
+                        STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
+                        dimensions.push_back(*dimensionsAsBitVector.back().begin());
+                        lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
+                        optimizationDirections.push_back(dirVar.first);
+                    }
+                    STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
+                    
+                    
+                    MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None)));
+                        
+                    // initialize data that will be needed for each epoch
+                    auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+                    auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+                    std::vector<ValueType> x, b;
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+                    
+                    if (currentEpochValues[0] < 0 && currentEpochValues[1] < 0) {
+                        // This case can only happen in these cases:
+                        assert(lowerCostBounds[0]);
+                        assert(currentEpochValues[0] == -1);
+                        assert(currentEpochValues[1] == -1);
+                        // This case has been checked already, so we can skip it.
+                        // Skipping this is actually necessary, since the rewardUnfolding will handle formulas like [F{"a"}>A,{"b"}>B "init"] incorrectly if A=B=-1.
+                        ++currentEpochValues[1];
+                    }
+                    std::set<typename EpochManager::Epoch> exploredEpochs;
+                    while (true) {
+                        auto currEpoch = rewardUnfolding.getStartEpoch(true);
+                        for (uint64_t i = 0; i < 2; ++i) {
+                            if (currentEpochValues[i] >= 0) {
+                                rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) currentEpochValues[i]);
+                            } else {
+                                rewardUnfolding.getEpochManager().setBottomDimension(currEpoch, dimensions[i]);
+                            }
+                        }
+                        auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
+                        for (auto const& epoch : epochOrder) {
+                            if (exploredEpochs.count(epoch) == 0) {
+                                auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
+                                rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
+                                ++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<ValueType> convertedPoint;
+                            for (uint64_t i = 0; i < 2; ++i) {
+                                if (lowerCostBounds[i]) {
+                                    // Translate > to >=
+                                    ++point[i];
+                                }
+                                if (i == 1 && storm::solver::maximize(optimizationDirections[i])) {
+                                    // When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction
+                                    --point[i];
+                                }
+                                if (point[i] < 0) {
+                                    // Skip this point
+                                    convertedPoint.clear();
+                                    continue;
+                                }
+                                convertedPoint.push_back(storm::utility::convertNumber<ValueType>(point[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor);
+                            }
+                            if (!convertedPoint.empty()) {
+                                resultPoints.push_back(std::move(convertedPoint));
+                            }
+                            
+                            if (currentEpochValues[1] == startEpochValues[1].first) {
+                                break;
+                            } else {
+                                ++currentEpochValues[0];
+                                currentEpochValues[1] = startEpochValues[1].first;
+                            }
+                        } else {
+                            ++currentEpochValues[1];
+                        }
+                    }
+                    
+    
+                    // When maximizing, there are border cases where one dimension can be arbitrarily large
+                    for (uint64_t i = 0; i < 2; ++i) {
+                        if (storm::solver::maximize(optimizationDirections[i]) && startEpochValues[i].first > 0) {
+                            std::vector<ValueType> newResultPoint(2);
+                            newResultPoint[i] = storm::utility::convertNumber<ValueType>(startEpochValues[i].first - 1) * startEpochValues[i].second;
+                            newResultPoint[1-i] = storm::utility::infinity<ValueType>();
+                            resultPoints.push_back(newResultPoint);
+                        }
+                    }
+                    
+                    filterDominatedPoints(resultPoints, optimizationDirections);
+                    return true;
+                }
+                
+                template<typename ModelType>
+                boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantileForDimension(Environment const& env, uint64_t dimension) const {
+                    // We assume that there is one bound value that violates the quantile and one bound value that satisfies it.
+                    
+                    // Let all other open bounds approach infinity
+                    std::set<storm::expressions::Variable> infinityVariables;
+                    for (auto const& d : getOpenDimensions()) {
+                        if (d != dimension) {
+                            infinityVariables.insert(getVariableForDimension(d));
+                        }
+                    }
+                    auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None));
+                    MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
+                    
+                    bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
+                    bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
+                    
+                    // initialize data that will be needed for each epoch
+                    auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+                    auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+                    std::vector<ValueType> x, b;
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
+                    
+                    uint64_t epochValue = 0;
+                    std::set<typename EpochManager::Epoch> 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<typename ModelType>
+                bool QuantileHelper<ModelType>::checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const {
+                    auto const& probabilityBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound();
+                    // Increase the precision until we get a conclusive result
+                    while (true) {
+                        ValueType numericResult = computeLimitValue(env, infDimensions);
+                        if (env.solver().isForceSoundness()) {
+                            auto lowerUpper = getLowerUpperBound(env, numericResult);
+                            bool lowerSat = probabilityBound.isSatisfied(lowerUpper.first);
+                            bool upperSat = probabilityBound.isSatisfied(lowerUpper.second);
+                            if (lowerSat == upperSat) {
+                                return lowerSat;
+                            } else {
+                                increasePrecision(env);
+                                ++numPrecisionRefinements;
+                            }
+                        } else {
+                            return probabilityBound.isSatisfied(numericResult);
+                        }
+                    }
+                }
+                
+                template<typename ModelType>
+                typename ModelType::ValueType QuantileHelper<ModelType>::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<BoundTransformation> bts(getDimension(), BoundTransformation::None);
+                    std::set<storm::expressions::Variable> 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<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
+                    if (!rewardUnfolding.getProb1Objectives().empty()) {
+                        assert(rewardUnfolding.getProb1Objectives().size() == 1);
+                        // The probability is one.
+                        return storm::utility::one<ValueType>();
+                    }
+                    // Get lower and upper bounds for the solution.
+                    auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
+                    auto upperBound = rewardUnfolding.getUpperObjectiveBound();
+
+                    // Initialize epoch models
+                    auto initEpoch = rewardUnfolding.getStartEpoch();
+                    auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
+                    // initialize data that will be needed for each epoch
+                    std::vector<ValueType> x, b;
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> 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<storm::models::sparse::Mdp<double>>;
+                template class QuantileHelper<storm::models::sparse::Mdp<storm::RationalNumber>>;
+
+            }
+        }
+    }
+}

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<typename ModelType>
+                class QuantileHelper {
+                    typedef typename ModelType::ValueType ValueType;
+                public:
+                    QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula);
+
+                    std::vector<std::vector<ValueType>> computeQuantile(Environment const& env) const;
+
+                private:
+                    std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env) const;
+                    bool exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const;
+
+                    /*!
+                     * Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
+                     */
+                    ValueType computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const;
+
+                    /*!
+                     * Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
+                     * The computed value is compared to the probability threshold.
+                     * In sound mode, precision is iteratively increased in case of 'inconsistent' results.
+                     */
+                    bool checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const;
+                    
+                    /// Computes the quantile with respect to the given dimension.
+                    /// boost::none is returned in case of insufficient precision.
+                    boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> computeQuantileForDimension(Environment const& env, uint64_t dim) const;
+                    
+                    /*!
+                     * Gets the number of dimensions of the underlying boudned until formula
+                     */
+                    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;
+                };
+            }
+        }
+    }
+}

src/storm/modelchecker/results/ExplicitParetoCurveCheckResult.h

@@ -12,8 +12,8 @@ namespace storm {
         class ExplicitParetoCurveCheckResult : public ParetoCurveCheckResult<ValueType> {
         public:
             ExplicitParetoCurveCheckResult();
-            ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector<typename ParetoCurveCheckResult<ValueType>::point_type> const& points, typename ParetoCurveCheckResult<ValueType>::polytope_type const& underApproximation, typename ParetoCurveCheckResult<ValueType>::polytope_type const& overApproximation);
-            ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector<typename ParetoCurveCheckResult<ValueType>::point_type>&& points, typename ParetoCurveCheckResult<ValueType>::polytope_type&& underApproximation, typename ParetoCurveCheckResult<ValueType>::polytope_type&& overApproximation);
+            ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector<typename ParetoCurveCheckResult<ValueType>::point_type> const& points, typename ParetoCurveCheckResult<ValueType>::polytope_type const& underApproximation = nullptr, typename ParetoCurveCheckResult<ValueType>::polytope_type const& overApproximation = nullptr);
+            ExplicitParetoCurveCheckResult(storm::storage::sparse::state_type const& state, std::vector<typename ParetoCurveCheckResult<ValueType>::point_type>&& points, typename ParetoCurveCheckResult<ValueType>::polytope_type&& underApproximation = nullptr, typename ParetoCurveCheckResult<ValueType>::polytope_type&& overApproximation = nullptr);
             
             ExplicitParetoCurveCheckResult(ExplicitParetoCurveCheckResult const& other) = default;
             ExplicitParetoCurveCheckResult& operator=(ExplicitParetoCurveCheckResult const& other) = default;

src/storm/modelchecker/results/ParetoCurveCheckResult.cpp

@@ -31,22 +31,38 @@ namespace storm {
             return points;
         }
         
+        template<typename ValueType>
+        bool ParetoCurveCheckResult<ValueType>::hasUnderApproximation() const {
+            return bool(underApproximation);
+        }
+        
+        template<typename ValueType>
+        bool ParetoCurveCheckResult<ValueType>::hasOverApproximation() const {
+            return bool(overApproximation);
+        }
+        
         template<typename ValueType>
         typename ParetoCurveCheckResult<ValueType>::polytope_type const& ParetoCurveCheckResult<ValueType>::getUnderApproximation() const {
+            STORM_LOG_ASSERT(hasUnderApproximation(), "Requested under approx. of Pareto curve although it does not exist.");
             return underApproximation;
         }
         
         template<typename ValueType>
         typename ParetoCurveCheckResult<ValueType>::polytope_type const& ParetoCurveCheckResult<ValueType>::getOverApproximation() const {
+            STORM_LOG_ASSERT(hasUnderApproximation(), "Requested over approx. of Pareto curve although it does not exist.");
             return overApproximation;
         }
         
         template<typename ValueType>
         std::ostream& ParetoCurveCheckResult<ValueType>::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){

src/storm/modelchecker/results/ParetoCurveCheckResult.h

@@ -19,14 +19,16 @@ namespace storm {
             virtual bool isParetoCurveCheckResult() const override;
             
             std::vector<point_type> 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<point_type> const& points, polytope_type const& underApproximation, polytope_type const& overApproximation);
-            ParetoCurveCheckResult(std::vector<point_type>&& points, polytope_type&& underApproximation, polytope_type&& overApproximation);
+            ParetoCurveCheckResult(std::vector<point_type> const& points, polytope_type const& underApproximation = nullptr, polytope_type const& overApproximation = nullptr);
+            ParetoCurveCheckResult(std::vector<point_type>&& points, polytope_type&& underApproximation = nullptr, polytope_type&& overApproximation = nullptr);
 
             // The pareto optimal points that have been found.
             std::vector<point_type> points;

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<storm::RationalFunctionVariable> getAllParameters(Model<storm::RationalFunction> const& model) {
+                std::set<storm::RationalFunctionVariable> parameters = getProbabilityParameters(model);
+                std::set<storm::RationalFunctionVariable> rewardParameters = getRewardParameters(model);
+                parameters.insert(rewardParameters.begin(), rewardParameters.end());
+                return parameters;
+            }
 #endif
-            
+
             template class Model<double>;
             template class Model<float>;
 

src/storm/models/sparse/Model.h

@@ -419,10 +419,28 @@ namespace storm {
                 boost::optional<std::shared_ptr<storm::storage::sparse::ChoiceOrigins>> choiceOrigins;
                 
             };
-            
+
 #ifdef STORM_HAVE_CARL
+            /*!
+             * Get all probability parameters occurring on transitions.
+             * @param model Model.
+             * @return Set of parameters.
+             */
             std::set<storm::RationalFunctionVariable> getProbabilityParameters(Model<storm::RationalFunction> const& model);
+
+            /*!
+             * Get all parameters occurring in rewards.
+             * @param model Model.
+             * @return Set of parameters.
+             */
             std::set<storm::RationalFunctionVariable> getRewardParameters(Model<storm::RationalFunction> const& model);
+
+            /*!
+             * Get all parameters (probability and rewards) occurring in the model.
+             * @param model Model.
+             * @return Set of parameters.
+             */
+            std::set<storm::RationalFunctionVariable> getAllParameters(Model<storm::RationalFunction> const& model);
 #endif
         } // namespace sparse
     } // namespace models

src/storm/storage/MaximalEndComponentDecomposition.cpp

@@ -74,10 +74,10 @@ namespace storm {
             if (states) {
                 endComponentStateSets.emplace_back(states->begin(), states->end(), true);
             } else {
-                std::vector<storm::storage::sparse::state_type> states;
-                states.resize(transitionMatrix.getRowGroupCount());
-                std::iota(states.begin(), states.end(), 0);
-                endComponentStateSets.emplace_back(states.begin(), states.end(), true);
+                std::vector<storm::storage::sparse::state_type> 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;

src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp

@@ -211,7 +211,7 @@ namespace storm {
             }
             
             // Finally construct the quotient model.
-            this->quotient = std::make_shared<ModelType>(builder.build(), std::move(newLabeling), std::move(rewardModels));
+            this->quotient = std::make_shared<ModelType>(builder.build(0,this->size(), this->size()), std::move(newLabeling), std::move(rewardModels));
         }
         
         template<typename ModelType>

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";

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;

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<ValueType> 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<ValueType> 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<uint_fast64_t> newRowGroupIndices;

src/storm/utility/graph.cpp

@@ -116,6 +116,43 @@ namespace storm {
                 return result;
             }
             
+            template <typename T>
+            bool hasCycle(storm::storage::SparseMatrix<T> const& transitionMatrix, boost::optional<storm::storage::BitVector> 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<uint64_t> 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 <typename T>
             bool checkIfECWithChoiceExists(storm::storage::SparseMatrix<T> const& transitionMatrix, storm::storage::SparseMatrix<T> 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<double> const& transitionMatrix);
            
+            template bool hasCycle(storm::storage::SparseMatrix<double> const& transitionMatrix, boost::optional<storm::storage::BitVector> const& subsystem);
+            
             template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices);
             
             template std::vector<uint_fast64_t> getDistances(storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional<storm::storage::BitVector> const& subsystem);
@@ -1727,8 +1766,10 @@ namespace storm {
             template storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix<storm::RationalNumber> 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<storm::storage::BitVector> const& choiceFilter);
             
             template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix);
-           
-           template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, storm::storage::SparseMatrix<storm::RationalNumber> const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices);
+            
+            template bool hasCycle(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, boost::optional<storm::storage::BitVector> const& subsystem);
+
+            template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, storm::storage::SparseMatrix<storm::RationalNumber> const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices);
 
             template std::vector<uint_fast64_t> getDistances(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional<storm::storage::BitVector> const& subsystem);
             
@@ -1785,6 +1826,8 @@ namespace storm {
             
             template storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix);
             
+            template bool hasCycle(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, boost::optional<storm::storage::BitVector> const& subsystem);
+
             template bool checkIfECWithChoiceExists(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, storm::storage::BitVector const& subsystem, storm::storage::BitVector const& choices);
             
             template std::vector<uint_fast64_t> getDistances(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, storm::storage::BitVector const& initialStates, boost::optional<storm::storage::BitVector> const& subsystem);

src/storm/utility/graph.h

@@ -81,6 +81,15 @@ namespace storm {
             template<typename T>
             storm::storage::BitVector getBsccCover(storm::storage::SparseMatrix<T> 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<typename T>
+            bool hasCycle(storm::storage::SparseMatrix<T> const& transitionMatrix, boost::optional<storm::storage::BitVector> const& subsystem = boost::none);
+            
             /*!
              * Checks whether there is an End Component that
              * 1. contains at least one of the specified choices and