Quantiles: made the SparseMdpPrctlModelChecker call the QuantileHelper for quantile formulas

src/storm/logic/Formula.cpp

@@ -209,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

@@ -110,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/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/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.computeMultiDimensionalQuantile();
+            
+            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/rewardbounded/Dimension.h

@@ -3,6 +3,7 @@
 #include <boost/optional.hpp>
 
 #include "storm/storage/BitVector.h"
+#include "storm/solver/OptimizationDirection.h"
 
 namespace storm {
     namespace modelchecker {
@@ -18,6 +19,7 @@ namespace storm {
                     ValueType scalingFactor;
                     storm::storage::BitVector dependentDimensions;
                     boost::optional<uint64_t> maxValue;
+                    boost::optional<storm::solver::OptimizationDirection> optimizationDirection;
                 };
             }
         }

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

@@ -88,6 +88,7 @@ namespace storm {
                                 dimension.memoryLabel = memLabel;
                                 dimension.isUpperBounded = subformula.hasUpperBound(dim);
                                 // for simplicity we do not allow intervals or unbounded formulas.
+                                // TODO: Quantiles: allow unbounded formulas
                                 STORM_LOG_THROW(subformula.hasLowerBound(dim) != dimension.isUpperBounded, storm::exceptions::NotSupportedException, "Bounded until formulas are only supported by this method if they consider either an upper bound or a lower bound. Got " << subformula << " instead.");
                                 // lower bounded until formulas with non-trivial left hand side are excluded as this would require some additional effort (in particular the ProductModel::transformMemoryState method).
                                 STORM_LOG_THROW(dimension.isUpperBounded || subformula.getLeftSubformula(dim).isTrueFormula(), storm::exceptions::NotSupportedException, "Lower bounded until formulas are only supported by this method if the left subformula is 'true'. Got " << subformula << " instead.");
@@ -216,20 +217,22 @@ namespace storm {
                             bound = dimFormula.asCumulativeRewardFormula().getBound();
                             isStrict = dimFormula.asCumulativeRewardFormula().isBoundStrict();
                         }
-                        STORM_LOG_THROW(!bound.containsVariables(), storm::exceptions::NotSupportedException, "The bound " << bound << " contains undefined constants.");
-                        ValueType discretizedBound = storm::utility::convertNumber<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()) {
+                            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>();
+                                }
+                            } else {
+                                discretizedBound = storm::utility::floor(discretizedBound);
                             }
-                        } 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;
                         }
-                        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;
                     }
                 }
         

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

@@ -0,0 +1,10 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+            }
+        }
+    }
+}

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

@@ -0,0 +1,21 @@
+#pragma once
+
+#include "storm/logic/QuantileFormula.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                template<typename ModelType>
+                class QuantileHelper {
+                    typedef typename ModelType::ValueType ValueType;
+                public:
+                    QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& formula) {}
+                    
+                    std::vector<std::vector<ValueType>> computeMultiDimensionalQuantile() { return {{27}};}
+                };
+            }
+        }
+    }
+}

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;