Added Mdp Region checking in storm.h, Some STORM_LOG_DEBUGs, fixes for sampling to work on Mdps

Former-commit-id: ab42fefd92

examples/pmdp/consensus/coin2_2.nm

@@ -0,0 +1,65 @@
+// COIN FLIPPING PROTOCOL FOR POLYNOMIAL RANDOMIZED CONSENSUS [AH90] 
+// gxn/dxp 20/11/00
+
+mdp
+
+// constants
+const int N=2;
+const int K=2;
+const int range = 2*(K+1)*N;
+const int counter_init = (K+1)*N;
+const int left = N;
+const int right = 2*(K+1)*N - N;
+
+//Coin Probabilities
+const double p;
+const double q;
+
+// shared coin
+global counter : [0..range] init counter_init;
+
+module process1
+	
+	// program counter
+	pc1 : [0..3];
+	// 0 - flip
+	// 1 - write 
+	// 2 - check
+	// 3 - finished
+	
+	// local coin
+	coin1 : [0..1];	
+
+	// flip coin
+	[] (pc1=0)  -> p : (coin1'=0) & (pc1'=1) + 1-p : (coin1'=1) & (pc1'=1);
+	// write tails -1  (reset coin to add regularity)
+	[] (pc1=1) & (coin1=0) & (counter>0) -> 1 : (counter'=counter-1) & (pc1'=2) & (coin1'=0);
+	// write heads +1 (reset coin to add regularity)
+	[] (pc1=1) & (coin1=1) & (counter<range) -> 1 : (counter'=counter+1) & (pc1'=2) & (coin1'=0);
+	// check
+	// decide tails
+	[] (pc1=2) & (counter<=left) -> 1 : (pc1'=3) & (coin1'=0);
+	// decide heads
+	[] (pc1=2) & (counter>=right) -> 1 : (pc1'=3) & (coin1'=1);
+	// flip again
+	[] (pc1=2) & (counter>left) & (counter<right) -> 1 : (pc1'=0);
+	// loop (all loop together when done)
+	[done] (pc1=3) -> 1 : (pc1'=3);
+
+endmodule
+
+// construct remaining processes through renaming
+module process2 = process1[pc1=pc2,coin1=coin2,p=q] endmodule
+
+// labels
+label "finished" = pc1=3 & pc2=3 ;
+label "all_coins_equal_0" = coin1=0 & coin2=0 ;
+label "all_coins_equal_1" = coin1=1 & coin2=1 ;
+label "agree" = coin1=coin2 ;
+label "finish_with_1" = pc1=3 & pc2=3 & coin1=1 & coin2=1 ;
+
+// rewards
+rewards "steps"
+	true : 1;
+endrewards
+

src/modelchecker/region/AbstractSparseRegionModelChecker.cpp

@@ -87,6 +87,7 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::specifyFormula(std::shared_ptr<storm::logic::Formula> formula) {
                 std::chrono::high_resolution_clock::time_point timeSpecifyFormulaStart = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Specifying the formula " << formula);
                 STORM_LOG_THROW(this->canHandle(*formula), storm::exceptions::InvalidArgumentException, "Tried to specify a formula that can not be handled.");
                 //Initialize the context for this formula
                 if (formula->isProbabilityOperatorFormula()) {
@@ -136,6 +137,7 @@ namespace storm {
                     }
                 } else if (this->isResultConstant() && this->constantResult.get() == storm::utility::region::convertNumber<ConstantType>(-1.0)){
                     //In this case, the result is constant but has not been computed yet. so do it now!
+                    STORM_LOG_DEBUG("The Result is constant and will be computed now.");
                     initializeSamplingModel(*this->getSimpleModel(), this->getSimpleFormula());
                     std::map<VariableType, CoefficientType> emptySubstitution;
                     this->getSamplingModel()->instantiate(emptySubstitution);
@@ -151,6 +153,7 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula) {
                 std::chrono::high_resolution_clock::time_point timeInitApproxModelStart = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("The Approximation Model is initialized");
                 STORM_LOG_THROW(this->isApproximationApplicable, storm::exceptions::UnexpectedException, "Approximation model requested but approximation is not applicable");
                 this->approximationModel=std::make_shared<ApproximationModel<ParametricSparseModelType, ConstantType>>(model, formula);
                 std::chrono::high_resolution_clock::time_point timeInitApproxModelEnd = std::chrono::high_resolution_clock::now();
@@ -160,6 +163,7 @@ namespace storm {
 
             template<typename ParametricSparseModelType, typename ConstantType>
             void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeSamplingModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula) {
+                STORM_LOG_DEBUG("The Sampling Model is initialized");
                 std::chrono::high_resolution_clock::time_point timeInitSamplingModelStart = std::chrono::high_resolution_clock::now();
                 this->samplingModel=std::make_shared<SamplingModel<ParametricSparseModelType, ConstantType>>(model, formula);
                 std::chrono::high_resolution_clock::time_point timeInitSamplingModelEnd = std::chrono::high_resolution_clock::now();
@@ -331,6 +335,7 @@ namespace storm {
 
             template<typename ParametricSparseModelType, typename ConstantType>
             void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::printStatisticsToStream(std::ostream& outstream) {
+                STORM_LOG_DEBUG("Printing statistics");
 
                 if(this->getSpecifiedFormula()==nullptr){
                     outstream << "Region Model Checker Statistics Error: No formula specified." << std::endl; 

src/modelchecker/region/ApproximationModel.h

@@ -45,7 +45,7 @@ namespace storm {
                 void instantiate(ParameterRegion<ParametricType> const& region);
 
                 /*!
-                 * Returns the approximated reachability probabilities for every state.
+                 * Returns the approximated reachability probabilities or expected values for every state.
                  * Undefined behavior if model has not been instantiated first!
                  * @param optimalityType Use MAXIMIZE to get upper bounds or MINIMIZE to get lower bounds
                  */

src/modelchecker/region/SamplingModel.cpp

@@ -98,13 +98,19 @@ namespace storm {
                                                                                                                     TableEntry* constantEntry) {
                 // Run through the rows of the original model and obtain the probability evaluation table, a matrix with dummy entries and the mapping between the two.
                 bool addRowGroups = parametricModel.getTransitionMatrix().hasNontrivialRowGrouping();
-                auto currRowGroup = parametricModel.getTransitionMatrix().getRowGroupIndices().begin();
-                storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getTransitionMatrix().getRowCount(), parametricModel.getTransitionMatrix().getColumnCount(), parametricModel.getTransitionMatrix().getEntryCount(), addRowGroups);
+                auto curRowGroup = parametricModel.getTransitionMatrix().getRowGroupIndices().begin();
+                storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getTransitionMatrix().getRowCount(),
+                                                                                parametricModel.getTransitionMatrix().getColumnCount(),
+                                                                                parametricModel.getTransitionMatrix().getEntryCount(),
+                                                                                true, //force dimensions
+                                                                                addRowGroups,
+                                                                                addRowGroups ? parametricModel.getTransitionMatrix().getRowGroupCount() : 0);
                 matrixEntryToEvalTableMapping.reserve(parametricModel.getTransitionMatrix().getEntryCount());
                 std::size_t numOfNonConstEntries=0;
                 for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
-                    if(addRowGroups && row==*currRowGroup){
+                    if(addRowGroups && row==*curRowGroup){
                         matrixBuilder.newRowGroup(row);
+                        ++curRowGroup;
                     }
                     ConstantType dummyEntry=storm::utility::one<ConstantType>();
                     for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
@@ -194,11 +200,12 @@ namespace storm {
                 }
                 std::unique_ptr<storm::modelchecker::CheckResult> resultPtr;
                 //perform model checking
+                boost::optional<storm::solver::OptimizationDirection> opDir = storm::logic::isLowerBound(this->formula->getComparisonType()) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize;
                 if(this->computeRewards){
                     resultPtr = modelChecker->computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula());
                 }
                 else {
-                    resultPtr = modelChecker->computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula());
+                    resultPtr = modelChecker->computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula(), false, opDir);
                 }
                 return resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
             }

src/modelchecker/region/SamplingModel.h

@@ -44,7 +44,7 @@ namespace storm {
                 void instantiate(std::map<VariableType, CoefficientType>const& point);
 
                 /*!
-                 * Returns the reachability probabilities for every state according to the current instantiation.
+                 * Returns the reachability probabilities (or the expected rewards) for every state according to the current instantiation.
                  * Undefined behavior if model has not been instantiated first!
                  */
                 std::vector<ConstantType> const& computeValues();

src/modelchecker/region/SparseDtmcRegionModelChecker.cpp

@@ -73,6 +73,7 @@ namespace storm {
                                                                                                    std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula,
                                                                                                    bool& isApproximationApplicable,
                                                                                                    boost::optional<ConstantType>& constantResult){
+                STORM_LOG_DEBUG("Preprocessing for DTMC started.");
                 STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidArgumentException, "Input model is required to have exactly one initial state.");
                 //Reset some data
                 this->smtSolver=nullptr;
@@ -88,6 +89,7 @@ namespace storm {
                 } else {
                     preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, constantResult);
                 }
+                STORM_LOG_DEBUG("Elimination of states with constant outgoing transitions is happening now.");
                 // Determine the set of states that is reachable from the initial state without jumping over a target state.
                 storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->getModel().getTransitionMatrix(), this->getModel().getInitialStates(), maybeStates, targetStates);
                 // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
@@ -145,6 +147,7 @@ namespace storm {
                 STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions.");
 
                 //Build the simple model
+                STORM_LOG_DEBUG("Building the resulting simplified model.");
                 //The matrix. The flexibleTransitions matrix might have empty rows where states have been eliminated.
                 //The new matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
                 std::vector<storm::storage::sparse::state_type> newStateIndexMap(flexibleTransitions.getNumberOfRows(), flexibleTransitions.getNumberOfRows()); //initialize with some illegal index to easily check if a transition leads to an unselected state
@@ -218,6 +221,7 @@ namespace storm {
                 // the final model
                 simpleModel = std::make_shared<ParametricSparseModelType>(matrixBuilder.build(), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
                 // the corresponding formula
+                STORM_LOG_DEBUG("Building the resulting simplified formula.");
                 std::shared_ptr<storm::logic::AtomicLabelFormula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
                 if(this->isComputeRewards()){
                     std::shared_ptr<storm::logic::ReachabilityRewardFormula> reachRewFormula(new storm::logic::ReachabilityRewardFormula(targetFormulaPtr));
@@ -238,6 +242,7 @@ namespace storm {
                                                                                                                    storm::storage::BitVector& targetStates,
                                                                                                                    bool& isApproximationApplicable,
                                                                                                                    boost::optional<ConstantType>& constantResult) {
+                STORM_LOG_DEBUG("Preprocessing for Dtmcs and reachability probabilities invoked.");
                 //Get Target States
                 storm::logic::AtomicLabelFormula const& labelFormula = this->getSpecifiedFormula()->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula().getSubformula().asAtomicLabelFormula();
                 storm::modelchecker::SparsePropositionalModelChecker<ParametricSparseModelType> modelChecker(this->getModel());
@@ -283,6 +288,7 @@ namespace storm {
                                                                                                              std::vector<ParametricType>& stateRewards,
                                                                                                              bool& isApproximationApplicable,
                                                                                                              boost::optional<ConstantType>& constantResult) {
+                STORM_LOG_DEBUG("Preprocessing for Dtmcs and reachability rewards invoked.");
                 //get the correct reward model
                 ParametricRewardModelType const* rewardModel;
                 if(this->getSpecifiedFormula()->asRewardOperatorFormula().hasRewardModelName()){
@@ -372,6 +378,7 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::computeReachabilityFunction(ParametricSparseModelType const& simpleModel){
                 std::chrono::high_resolution_clock::time_point timeComputeReachabilityFunctionStart = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Computing the Reachability function...");
                 //get the one step probabilities and the transition matrix of the simple model without target/sink state
                 storm::storage::SparseMatrix<ParametricType> backwardTransitions(simpleModel.getBackwardTransitions());
                 std::vector<ParametricType> oneStepProbabilities(simpleModel.getNumberOfStates()-2, storm::utility::zero<ParametricType>());
@@ -397,7 +404,7 @@ namespace storm {
                                     this->reachabilityFunction->denominator().toString(false, true) +
                                 " )";
                     std::cout << std::endl <<"the resulting reach prob function is " << std::endl << funcStr << std::endl << std::endl;*/
-                STORM_LOG_DEBUG("Computed reachabilityFunction");
+                STORM_LOG_DEBUG("Done with computing the reachabilityFunction");
                 std::chrono::high_resolution_clock::time_point timeComputeReachabilityFunctionEnd = std::chrono::high_resolution_clock::now();
                 this->timeComputeReachabilityFunction = timeComputeReachabilityFunctionEnd-timeComputeReachabilityFunctionStart;
             }
@@ -653,6 +660,7 @@ namespace storm {
 
             template<typename ParametricSparseModelType, typename ConstantType>
             void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeSMTSolver() {
+                STORM_LOG_DEBUG("Initializing the Smt Solver");
 
                 storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions..
                 this->smtSolver = std::shared_ptr<storm::solver::Smt2SmtSolver>(new storm::solver::Smt2SmtSolver(manager, true));

src/modelchecker/region/SparseMdpRegionModelChecker.cpp

@@ -72,16 +72,22 @@ namespace storm {
                                                                                                   std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula,
                                                                                                   bool& isApproximationApplicable,
                                                                                                   boost::optional<ConstantType>& constantResult){
+                STORM_LOG_DEBUG("Preprocessing for MDPs started.");
                 STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidArgumentException, "Input model is required to have exactly one initial state.");
                 storm::storage::BitVector maybeStates, targetStates;
                 preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, constantResult);
                 //TODO: Actually get a more simple model. This makes a deep copy of the model...
-                simpleModel = std::make_shared<ParametricSparseModelType>(this->getModel());
+                STORM_LOG_WARN("No simplification of the original model (like elimination of constant transitions) is happening. Will just use a copy of the original model");
+                simpleModel = std::make_shared<ParametricSparseModelType>(this->getModel()); //Note: an actual copy is technically not necessary.. but we will do it here..
                 simpleFormula = this->getSpecifiedFormula();
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            void SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult) {
+            void SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForProbabilities(storm::storage::BitVector& maybeStates,
+                                                                                                                  storm::storage::BitVector& targetStates,
+                                                                                                                  bool& isApproximationApplicable,
+                                                                                                                  boost::optional<ConstantType>& constantResult) {
+                STORM_LOG_DEBUG("Preprocessing for Mdps and reachability probabilities invoked.");
                 //Get Target States
                 storm::logic::AtomicLabelFormula const& labelFormula = this->getSpecifiedFormula()->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula().getSubformula().asAtomicLabelFormula();
                 storm::modelchecker::SparsePropositionalModelChecker<ParametricSparseModelType> modelChecker(this->getModel());
@@ -246,6 +252,11 @@ namespace storm {
                 return false;
             }
             
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::checkSmt(ParameterRegion<ParametricType>& region) {
+                STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "checkSmt invoked but smt solving has not been implemented for MDPs.");
+            }
+
 #ifdef STORM_HAVE_CARL
             template class SparseMdpRegionModelChecker<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
 #endif

src/modelchecker/region/SparseMdpRegionModelChecker.h

@@ -80,6 +80,16 @@ namespace storm {
                  */
                 virtual bool checkPoint(ParameterRegion<ParametricType>& region, std::map<VariableType, CoefficientType>const& point, bool favorViaFunction=false);
                 
+                /*!
+                 * Starts the SMTSolver to get the result.
+                 * The current regioncheckresult of the region should be EXISTSSAT or EXISTVIOLATED.
+                 * Otherwise, a sampingPoint will be computed.
+                 * True is returned iff the solver was successful (i.e., it returned sat or unsat)
+                 * A Sat- or Violated point is set, if the solver has found one (not yet implemented!).
+                 * The region checkResult of the given region is changed accordingly.
+                 */
+                virtual bool checkSmt(ParameterRegion<ParametricType>& region); 
+
             };
         } //namespace region
     } // namespace modelchecker

src/utility/region.h

@@ -20,14 +20,6 @@
 #include "src/logic/ComparisonType.h"
 #include <src/adapters/CarlAdapter.h>
 
-// Forward-declare region modelchecker  class.
-    namespace storm {
-        namespace modelchecker{
-        template<typename ParametricType, typename ConstantType>
-        class SparseDtmcRegionModelChecker;
-        }
-    }
-
 namespace storm {
     namespace utility{
         namespace region {

src/utility/storm.h

@@ -54,6 +54,7 @@
 #include "src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h"
 #include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
 #include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
+#include "src/modelchecker/region/SparseMdpRegionModelChecker.h"
 #include "src/modelchecker/region/ParameterRegion.h"
 #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
@@ -243,27 +244,32 @@ namespace storm {
     inline void verifySparseModel(std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
 
         for (auto const& formula : formulas) {
-            STORM_LOG_THROW(model->getType() == storm::models::ModelType::Dtmc, storm::exceptions::InvalidSettingsException, "Currently parametric verification is only available for DTMCs.");
-            std::shared_ptr<storm::models::sparse::Dtmc<storm::RationalFunction>> dtmc = model->template as<storm::models::sparse::Dtmc<storm::RationalFunction>>();
-            
             if(storm::settings::generalSettings().isParametricRegionSet()){
                 std::cout << std::endl << "Model checking property: " << *formula << " for all parameters in the given regions." << std::endl;
-
                 auto regions=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::getRegionsFromSettings();
+                if(model->getType() == storm::models::ModelType::Dtmc){
+                    std::shared_ptr<storm::models::sparse::Dtmc<storm::RationalFunction>> dtmc = model->template as<storm::models::sparse::Dtmc<storm::RationalFunction>>();
                     storm::modelchecker::region::SparseDtmcRegionModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, double> modelchecker(*dtmc);
-                if (modelchecker.canHandle(*formula.get())) {
+                    STORM_LOG_THROW(modelchecker.canHandle(*formula.get()), storm::exceptions::InvalidSettingsException, "The parametric region check engine currently does not support this property.");
                     modelchecker.specifyFormula(formula);
                     modelchecker.checkRegions(regions);
-                }
-                else {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The parametric region check engine currently does not support this property.");
+                    modelchecker.printStatisticsToStream(std::cout);
+                } else if (model->getType() == storm::models::ModelType::Mdp){
+                    std::shared_ptr<storm::models::sparse::Mdp<storm::RationalFunction>> mdp = model->template as<storm::models::sparse::Mdp<storm::RationalFunction>>();
+                    storm::modelchecker::region::SparseMdpRegionModelChecker<storm::models::sparse::Mdp<storm::RationalFunction>, double> modelchecker(*mdp);
+                    STORM_LOG_THROW(modelchecker.canHandle(*formula.get()), storm::exceptions::InvalidSettingsException, "The parametric region check engine currently does not support this property.");
+                    modelchecker.specifyFormula(formula);
+                    modelchecker.checkRegions(regions);
+                    modelchecker.printStatisticsToStream(std::cout);
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Currently parametric region verification is only available for DTMCs and Mdps.");
                 }
                // for(auto const& reg : regions){
                //     std::cout << reg.toString() << "      Result: " << reg.getCheckResult() << std::endl;
                // }
-                modelchecker.printStatisticsToStream(std::cout);
-
             } else {
+                STORM_LOG_THROW(model->getType() == storm::models::ModelType::Dtmc, storm::exceptions::InvalidSettingsException, "Currently parametric verification via state elimination is only available for DTMCs.");
+                std::shared_ptr<storm::models::sparse::Dtmc<storm::RationalFunction>> dtmc = model->template as<storm::models::sparse::Dtmc<storm::RationalFunction>>();
                 std::cout << std::endl << "Model checking property: " << *formula << " ...";
                 std::unique_ptr<storm::modelchecker::CheckResult> result;
 
@@ -420,4 +426,3 @@ namespace storm {
 }
 
 #endif	/* STORM_H */
-