storm-pars: Added possibility to compute the extremal value within a given region using parameter lifting.

CHANGELOG.md

@@ -21,6 +21,7 @@ Version 1.2.x
 - `storm-gspn`: Added support for single/infinite/k-server semantics for GSPNs given in the .pnpro format
 - `storm-gspn`: Added option to set a global capacity for all places
 - `storm-gspn`: Added option to include a set of standard properties when converting GSPNs to jani
+- `storm-pars`: Added possibility to compute the extremal value within a given region using parameter lifting
 
 ### Version 1.2.3 (2018/07)
 - Fix in version parsing

src/storm-pars-cli/storm-pars.cpp

@@ -365,6 +365,35 @@ namespace storm {
             }
         }
         
+        template <typename ValueType>
+        void computeRegionExtremumWithSparseEngine(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, SymbolicInput const& input, std::vector<storm::storage::ParameterRegion<ValueType>> const& regions) {
+            STORM_LOG_ASSERT(!regions.empty(), "Can not analyze an empty set of regions.");
+            auto regionSettings = storm::settings::getModule<storm::settings::modules::RegionSettings>();
+            auto engine = regionSettings.getRegionCheckEngine();
+            storm::solver::OptimizationDirection direction = regionSettings.getExtremumDirection();
+            ValueType precision = storm::utility::convertNumber<ValueType>(regionSettings.getExtremumValuePrecision());
+            for (auto const& property : input.properties) {
+                for (auto const& region : regions) {
+                    STORM_PRINT_AND_LOG("Computing extremal value for property " << property.getName() << ": " << *property.getRawFormula() << " within region " << region << "..." << std::endl);
+                    storm::utility::Stopwatch watch(true);
+                    auto valueValuation = storm::api::computeExtremalValue<ValueType>(model, storm::api::createTask<ValueType>(property.getRawFormula(), true), region, engine, direction, precision);
+                    watch.stop();
+                    std::stringstream valuationStr;
+                    bool first = true;
+                    for (auto const& v : valueValuation.second) {
+                        if (first) {
+                            first = false;
+                        } else {
+                            valuationStr << ", ";
+                        }
+                        valuationStr << v.first << "=" << v.second;
+                    }
+                    STORM_PRINT_AND_LOG("Result at initial state: " << valueValuation.first << " ( approx. " << storm::utility::convertNumber<double>(valueValuation.first) << ") at [" << valuationStr.str() << "]." << std::endl)
+                    STORM_PRINT_AND_LOG("Time for model checking: " << watch << "." << std::endl);
+                }
+            }
+        }
+        
         template <typename ValueType>
         void verifyRegionsWithSparseEngine(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, SymbolicInput const& input, std::vector<storm::storage::ParameterRegion<ValueType>> const& regions) {
             STORM_LOG_ASSERT(!regions.empty(), "Can not analyze an empty set of regions.");
@@ -424,7 +453,12 @@ namespace storm {
             if (regions.empty()) {
                 storm::pars::verifyPropertiesWithSparseEngine(model, input, samples);
             } else {
-                storm::pars::verifyRegionsWithSparseEngine(model, input, regions);
+                auto regionSettings = storm::settings::getModule<storm::settings::modules::RegionSettings>();
+                if (regionSettings.isExtremumSet()) {
+                    storm::pars::computeRegionExtremumWithSparseEngine(model, input, regions);
+                } else {
+                    storm::pars::verifyRegionsWithSparseEngine(model, input, regions);
+                }
             }
         }
         

src/storm-pars/api/region.h

@@ -187,6 +187,20 @@ namespace storm {
             auto regionChecker = initializeRegionModelChecker(env, model, task, engine);
             return regionChecker->performRegionRefinement(env, region, coverageThreshold, refinementDepthThreshold, hypothesis);
         }
+    
+        /*!
+         * Finds the extremal value in the given region
+         * @param engine The considered region checking engine
+         * @param coverageThreshold if given, the refinement stops as soon as the fraction of the area of the subregions with inconclusive result is less then this threshold
+         * @param refinementDepthThreshold if given, the refinement stops at the given depth. depth=0 means no refinement.
+         * @param hypothesis if not 'unknown', it is only checked whether the hypothesis holds (and NOT the complementary result).
+         */
+        template <typename ValueType>
+        std::pair<ValueType, typename storm::storage::ParameterRegion<ValueType>::Valuation> computeExtremalValue(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, storm::modelchecker::CheckTask<storm::logic::Formula, ValueType> const& task, storm::storage::ParameterRegion<ValueType> const& region, storm::modelchecker::RegionCheckEngine engine, storm::solver::OptimizationDirection const& dir, boost::optional<ValueType> const& precision) {
+            Environment env;
+            auto regionChecker = initializeRegionModelChecker(env, model, task, engine);
+            return regionChecker->computeExtremalValue(env, region, dir, precision.is_initialized() ? precision.get() : storm::utility::zero<ValueType>());
+        }
         
         template <typename ValueType>
         void exportRegionCheckResultToFile(std::unique_ptr<storm::modelchecker::CheckResult> const& checkResult, std::string const& filename, bool onlyConclusiveResults = false) {

src/storm-pars/modelchecker/region/RegionModelChecker.cpp

@@ -152,6 +152,13 @@ namespace storm {
                 return std::make_unique<storm::modelchecker::RegionRefinementCheckResult<ParametricType>>(std::move(result), std::move(regionCopyForResult));
             }
 
+        template <typename ParametricType>
+        std::pair<ParametricType, typename storm::storage::ParameterRegion<ParametricType>::Valuation> RegionModelChecker<ParametricType>::computeExtremalValue(Environment const& env, storm::storage::ParameterRegion<ParametricType> const& region, storm::solver::OptimizationDirection const& dir, ParametricType const& precision) {
+            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing extremal values is not supported for this region model checker.");
+            return std::pair<ParametricType, typename storm::storage::ParameterRegion<ParametricType>::Valuation>();
+        }
+
+        
         template <typename ParametricType>
         bool RegionModelChecker<ParametricType>::isRegionSplitEstimateSupported() const {
             return false;

src/storm-pars/modelchecker/region/RegionModelChecker.h

@@ -59,6 +59,13 @@ namespace storm {
              */
             std::unique_ptr<storm::modelchecker::RegionRefinementCheckResult<ParametricType>> performRegionRefinement(Environment const& env, storm::storage::ParameterRegion<ParametricType> const& region, boost::optional<ParametricType> const& coverageThreshold, boost::optional<uint64_t> depthThreshold = boost::none, RegionResultHypothesis const& hypothesis = RegionResultHypothesis::Unknown);
             
+            /*!
+             * Finds the extremal value within the given region and with the given precision.
+             * The returned value v corresponds to the value at the returned valuation.
+             * The actual maximum (minimum) lies in the interval [v, v+precision] ([v-precision, v])
+             */
+            virtual std::pair<ParametricType, typename storm::storage::ParameterRegion<ParametricType>::Valuation> computeExtremalValue(Environment const& env, storm::storage::ParameterRegion<ParametricType> const& region, storm::solver::OptimizationDirection const& dir, ParametricType const& precision);
+            
             /*!
              * Returns true if region split estimation (a) was enabled when model and check task have been specified and (b) is supported by this region model checker.
              */

src/storm-pars/modelchecker/region/SparseParameterLiftingModelChecker.cpp

@@ -1,5 +1,7 @@
 #include "storm-pars/modelchecker/region/SparseParameterLiftingModelChecker.h"
 
+#include <queue>
+
 #include "storm/adapters/RationalFunctionAdapter.h"
 #include "storm/logic/FragmentSpecification.h"
 #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -143,6 +145,74 @@ namespace storm {
             return storm::utility::convertNumber<typename SparseModelType::ValueType>(getBound(env, region, dirForParameters)->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()]);
         }
         
+        template <typename SparseModelType, typename ConstantType>
+        struct RegionBound {
+            RegionBound(RegionBound<SparseModelType, ConstantType> const& other) = default;
+            RegionBound(storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& r, ConstantType const& b) : region(r), bound(b) {}
+            storm::storage::ParameterRegion<typename SparseModelType::ValueType> region;
+            ConstantType bound;
+        };
+        
+        template <typename SparseModelType, typename ConstantType>
+        std::pair<typename SparseModelType::ValueType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation> SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::computeExtremalValue(Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dir, typename SparseModelType::ValueType const& precision) {
+            STORM_LOG_THROW(this->parametricModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "Getting extremal values at the initial state requires a model with a single initial state.");
+            boost::optional<ConstantType> value;
+            typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation valuation;
+
+            for (auto const& v : region.getVerticesOfRegion(region.getVariables())) {
+                auto currValue = getInstantiationChecker().check(env, v)->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
+                if (!value.is_initialized() || (storm::solver::minimize(dir) ? currValue < value.get() : currValue > value.get())) {
+                    value = currValue;
+                    valuation = v;
+                    STORM_LOG_INFO("Current value for extremum: " << value.get() << ".");
+                }
+            }
+            
+            auto cmp = storm::solver::minimize(dir) ?
+                    [](RegionBound<SparseModelType, ConstantType> const& lhs, RegionBound<SparseModelType, ConstantType> const& rhs) { return lhs.bound > rhs.bound; } :
+                    [](RegionBound<SparseModelType, ConstantType> const& lhs, RegionBound<SparseModelType, ConstantType> const& rhs) { return lhs.bound < rhs.bound; };
+            std::priority_queue<RegionBound<SparseModelType, ConstantType>, std::vector<RegionBound<SparseModelType, ConstantType>>, decltype(cmp)> regionQueue(cmp);
+            regionQueue.emplace(region, storm::utility::zero<ConstantType>());
+            auto totalArea = storm::utility::convertNumber<ConstantType>(region.area());
+            auto coveredArea = storm::utility::zero<ConstantType>();
+            while (!regionQueue.empty()) {
+                auto const& currRegion = regionQueue.top().region;
+                
+                // Check whether this region contains a new 'good' value
+                auto currValue = getInstantiationChecker().check(env, currRegion.getCenterPoint())->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
+                if (storm::solver::minimize(dir) ? currValue < value.get() : currValue > value.get()) {
+                    value = currValue;
+                    valuation = currRegion.getCenterPoint();
+                }
+
+                // Check whether this region needs further investigation (i.e. splitting)
+                auto currBound = getBound(env, currRegion, dir)->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
+                std::vector<storm::storage::ParameterRegion<typename SparseModelType::ValueType>> newRegions;
+                if (storm::solver::minimize(dir)) {
+                    if (currBound < value.get() - storm::utility::convertNumber<ConstantType>(precision)) {
+                        currRegion.split(currRegion.getCenterPoint(), newRegions);
+                    }
+                } else {
+                    if (currBound > value.get() + storm::utility::convertNumber<ConstantType>(precision)) {
+                        currRegion.split(currRegion.getCenterPoint(), newRegions);
+                    }
+                }
+                
+                if (newRegions.empty()) {
+                    coveredArea += storm::utility::convertNumber<ConstantType>(currRegion.area());
+                }
+                regionQueue.pop();
+                for (auto const& r : newRegions) {
+                    regionQueue.emplace(r, currBound);
+                }
+                STORM_LOG_INFO("Current value : " << value.get() << ", current bound: " << regionQueue.top().bound << ".");
+                STORM_LOG_INFO("Covered " << (coveredArea * storm::utility::convertNumber<ConstantType>(100.0) / totalArea) << "% of the region.");
+            }
+            
+            return std::make_pair(storm::utility::convertNumber<typename SparseModelType::ValueType>(value.get()), valuation);
+        }
+
+        
         template <typename SparseModelType, typename ConstantType>
         SparseModelType const& SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getConsideredParametricModel() const {
             return *parametricModel;

src/storm-pars/modelchecker/region/SparseParameterLiftingModelChecker.h

@@ -49,6 +49,13 @@ namespace storm {
             virtual typename SparseModelType::ValueType getBoundAtInitState(Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters) override;
 
             
+            /*!
+             * Finds the extremal value within the given region and with the given precision.
+             * The returned value v corresponds to the value at the returned valuation.
+             * The actual maximum (minimum) lies in the interval [v, v+precision] ([v-precision, v])
+             */
+            virtual std::pair<typename SparseModelType::ValueType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation> computeExtremalValue(Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters, typename SparseModelType::ValueType const& precision) override;
+            
             SparseModelType const& getConsideredParametricModel() const;
             CheckTask<storm::logic::Formula, ConstantType> const& getCurrentCheckTask() const;
             

src/storm-pars/settings/modules/RegionSettings.cpp

@@ -20,6 +20,7 @@ namespace storm {
             const std::string RegionSettings::hypothesisOptionName = "hypothesis";
             const std::string RegionSettings::hypothesisShortOptionName = "hyp";
             const std::string RegionSettings::refineOptionName = "refine";
+            const std::string RegionSettings::extremumOptionName = "extremum";
             const std::string RegionSettings::checkEngineOptionName = "engine";
             const std::string RegionSettings::printNoIllustrationOptionName = "noillustration";
             const std::string RegionSettings::printFullResultOptionName = "printfullresult";
@@ -36,6 +37,11 @@ namespace storm {
                                 .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("coverage-threshold", "Refinement converges if the fraction of unknown area falls below this threshold.").setDefaultValueDouble(0.05).addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0,1.0)).build())
                                 .addArgument(storm::settings::ArgumentBuilder::createIntegerArgument("depth-limit", "If given, limits the number of times a region is refined.").setDefaultValueInteger(-1).build()).build());
                 
+                std::vector<std::string> directions = {"min", "max"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, extremumOptionName, false, "Computes the extremum within the region.")
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("direction", "The optimization direction").addValidatorString(storm::settings::ArgumentValidatorFactory::createMultipleChoiceValidator(directions)).build())
+                                .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("precision", "The desired precision").setDefaultValueDouble(0.05).addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0,1.0)).build()).build());
+                
                 std::vector<std::string> engines = {"pl", "exactpl", "validatingpl"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, checkEngineOptionName, true, "Sets which engine is used for analyzing regions.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the engine to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(engines)).setDefaultValueString("pl").build()).build());
@@ -93,6 +99,24 @@ namespace storm {
                 return (uint64_t) depth;
             }
             
+            bool RegionSettings::isExtremumSet() const {
+                return this->getOption(extremumOptionName).getHasOptionBeenSet();
+            }
+				
+            storm::solver::OptimizationDirection RegionSettings::getExtremumDirection() const {
+                auto str = this->getOption(extremumOptionName).getArgumentByName("direction").getValueAsString();
+                if (str == "min") {
+                    return storm::solver::OptimizationDirection::Minimize;
+                } else {
+                    assert(str == "max");
+                    return storm::solver::OptimizationDirection::Maximize;
+                }
+            }
+				
+            double RegionSettings::getExtremumValuePrecision() const {
+                return this->getOption(extremumOptionName).getArgumentByName("precision").getValueAsDouble();
+            }
+
             storm::modelchecker::RegionCheckEngine RegionSettings::getRegionCheckEngine() const {
                 std::string engineString = this->getOption(checkEngineOptionName).getArgumentByName("name").getValueAsString();
                 
@@ -110,6 +134,12 @@ namespace storm {
                 return result;
             }
             
+            bool RegionSettings::check() const {
+                if (isRefineSet() && isExtremumSet()) {
+                    STORM_LOG_ERROR("Can not compute extremum values AND perform region refinement.");
+                }
+            }
+
             bool RegionSettings::isPrintNoIllustrationSet() const {
                 return this->getOption(printNoIllustrationOptionName).getHasOptionBeenSet();
             }

src/storm-pars/settings/modules/RegionSettings.h

@@ -2,6 +2,7 @@
 
 #include "storm-pars/modelchecker/region/RegionCheckEngine.h"
 #include "storm-pars/modelchecker/region/RegionResultHypothesis.h"
+#include "storm/solver/OptimizationDirection.h"
 
 #include "storm/settings/modules/ModuleSettings.h"
 
@@ -61,6 +62,21 @@ namespace storm {
                  */
                 uint64_t getDepthLimit() const;
                 
+                /*!
+				 * Retrieves whether an extremal value is to be computed
+				 */
+				bool isExtremumSet() const;
+				
+				/*!
+				 * Retrieves whether to minimize or maximize when computing the extremal value
+				 */
+				storm::solver::OptimizationDirection getExtremumDirection() const;
+				
+				/*!
+				 * Retrieves the precision for the extremal value
+				 */
+				double getExtremumValuePrecision() const;
+
 				/*!
 				 * Retrieves which type of region check should be performed
 				 */
@@ -76,6 +92,8 @@ namespace storm {
                  */
                 bool isPrintFullResultSet() const;
                 
+                bool check() const override;
+                
                 const static std::string moduleName;
                 
             private:
@@ -84,6 +102,7 @@ namespace storm {
 				const static std::string hypothesisOptionName;
 				const static std::string hypothesisShortOptionName;
 				const static std::string refineOptionName;
+				const static std::string extremumOptionName;
 				const static std::string checkEngineOptionName;
 				const static std::string printNoIllustrationOptionName;
 				const static std::string printFullResultOptionName;

src/storm-pars/storage/ParameterRegion.h

@@ -16,8 +16,9 @@ namespace storm {
             ParameterRegion();
             ParameterRegion(Valuation const& lowerBoundaries, Valuation const& upperBoundaries);
             ParameterRegion(Valuation&& lowerBoundaries, Valuation&& upperBoundaries);
-            ParameterRegion(ParameterRegion const& other) = default;
-            ParameterRegion(ParameterRegion&& other) = default;
+            ParameterRegion(ParameterRegion<ParametricType> const& other) = default;
+            ParameterRegion(ParameterRegion<ParametricType>&& other) = default;
+            ParameterRegion<ParametricType>& operator=(ParameterRegion<ParametricType> const& other) = default;
             
             virtual ~ParameterRegion() = default;