PLA for continuous models

src/storm/modelchecker/parametric/ParameterLifting.cpp

@@ -81,7 +81,8 @@ namespace storm {
                 parameterLiftingCheckerStopwatch.start();
                 if(result == RegionCheckResult::ExistsSat || result == RegionCheckResult::CenterSat) {
                     // show AllSat:
-                    if(parameterLiftingChecker->check(region, this->currentCheckTask->getOptimizationDirection())->asExplicitQualitativeCheckResult()[*getConsideredParametricModel().getInitialStates().begin()]) {
+                    storm::solver::OptimizationDirection parameterOptimizationDirection = isLowerBound(this->currentCheckTask->getBound().comparisonType) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize;
+                    if(parameterLiftingChecker->check(region, parameterOptimizationDirection)->asExplicitQualitativeCheckResult()[*getConsideredParametricModel().getInitialStates().begin()]) {
                         result = RegionCheckResult::AllSat;
                     } else if (sampleVerticesOfRegion) {
                         parameterLiftingCheckerStopwatch.stop(); instantiationCheckerStopwatch.start();
@@ -96,7 +97,8 @@ namespace storm {
                     }
                 } else if (result == RegionCheckResult::ExistsViolated || result == RegionCheckResult::CenterViolated) {
                     // show AllViolated:
-                    if(!parameterLiftingChecker->check(region, storm::solver::invert(this->currentCheckTask->getOptimizationDirection()))->asExplicitQualitativeCheckResult()[*getConsideredParametricModel().getInitialStates().begin()]) {
+                    storm::solver::OptimizationDirection parameterOptimizationDirection = isLowerBound(this->currentCheckTask->getBound().comparisonType) ? storm::solver::OptimizationDirection::Maximize : storm::solver::OptimizationDirection::Minimize;
+                    if(!parameterLiftingChecker->check(region, parameterOptimizationDirection)->asExplicitQualitativeCheckResult()[*getConsideredParametricModel().getInitialStates().begin()]) {
                         result = RegionCheckResult::AllViolated;
                     } else if (sampleVerticesOfRegion) {
                         parameterLiftingCheckerStopwatch.stop(); instantiationCheckerStopwatch.start();

src/storm/models/sparse/Ctmc.cpp

@@ -36,6 +36,11 @@ namespace storm {
                 return exitRates;
             }
             
+            template <typename ValueType, typename RewardModelType>
+            std::vector<ValueType>& Ctmc<ValueType, RewardModelType>::getExitRateVector() {
+                return exitRates;
+            }
+            
             template <typename ValueType, typename RewardModelType>
             std::vector<ValueType> Ctmc<ValueType, RewardModelType>::createExitRateVector(storm::storage::SparseMatrix<ValueType> const& rateMatrix) {
                 std::vector<ValueType> exitRates(rateMatrix.getRowCount());

src/storm/models/sparse/Ctmc.h

@@ -65,6 +65,7 @@ namespace storm {
                  * @return The exit rate vector.
                  */
                 std::vector<ValueType> const& getExitRateVector() const;
+                std::vector<ValueType>& getExitRateVector();
                 
             private:
                 /*!

src/storm/storage/ParameterRegion.cpp

@@ -188,7 +188,7 @@ namespace storm {
                     var = std::make_unique<VariableType>(v);
                 }
             }
-            STORM_LOG_ASSERT(var, "Could not find parameter " << parameter << "in the set of considered variables");
+            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));

src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp

@@ -0,0 +1,149 @@
+#include "storm/transformer/ContinuousToDiscreteTimeModelTransformer.h"
+
+#include <unordered_map>
+
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/logic/Formulas.h"
+#include "storm/logic/FragmentSpecification.h"
+#include "storm/logic/CloneVisitor.h"
+#include "storm/utility/macros.h"
+#include "storm/utility/vector.h"
+
+#include "storm/exceptions/InvalidArgumentException.h"
+
+namespace storm {
+    namespace transformer {
+        
+        template <typename ValueType>
+        void transformContinuousToDiscreteModelInPlace(std::shared_ptr<storm::models::sparse::Model<ValueType>>& markovModel, std::shared_ptr<storm::logic::Formula const>& formula) {
+            boost::optional<std::string> timeRewardModelName;
+            if (formula->isTimeOperatorFormula()) {
+                auto const& timeOpFormula = formula->asTimeOperatorFormula();
+                if (timeOpFormula.getSubformula().isReachabilityTimeFormula()) {
+                    auto reachabilityRewardFormula = std::make_shared<storm::logic::EventuallyFormula>(storm::logic::CloneVisitor().clone(timeOpFormula.getSubformula().asReachabilityTimeFormula().getSubformula()), storm::logic::FormulaContext::Reward);
+                    timeRewardModelName = "time";
+                    // make sure that the reward model name is not already in use
+                    while (markovModel->hasRewardModel(*timeRewardModelName)) *timeRewardModelName += "_";
+                    formula = std::make_shared<storm::logic::RewardOperatorFormula const>(reachabilityRewardFormula, timeRewardModelName, timeOpFormula.getOperatorInformation());
+                }
+            }
+        
+            if (markovModel->isOfType(storm::models::ModelType::Ctmc)) {
+                SparseCtmcToSparseDtmcTransformer<storm::models::sparse::Ctmc<ValueType>> transformer;
+                if (transformer.transformationPreservesProperty(*formula)) {
+                    STORM_LOG_INFO("Transforming Ctmc to embedded Dtmc...");
+                    markovModel = transformer.translate(std::move(*markovModel->template as<storm::models::sparse::Ctmc<ValueType>>()), timeRewardModelName);
+                }
+            } else if (markovModel->isOfType(storm::models::ModelType::MarkovAutomaton)) {
+                SparseMaToSparseMdpTransformer<storm::models::sparse::MarkovAutomaton<ValueType>> transformer;
+                if (transformer.transformationPreservesProperty(*formula)) {
+                    STORM_LOG_INFO("Transforming Markov automaton to embedded Mdp...");
+                    markovModel = transformer.translate(std::move(*markovModel->template as<storm::models::sparse::MarkovAutomaton<ValueType>>()), timeRewardModelName);
+                }
+            }
+        }
+        
+        template<typename CtmcType>
+        std::shared_ptr<storm::models::sparse::Dtmc<typename CtmcType::ValueType, typename CtmcType::RewardModelType>> SparseCtmcToSparseDtmcTransformer<CtmcType>::translate(CtmcType&& ctmc, boost::optional<std::string> const& timeRewardModelName) {
+            // Turn the rates into probabilities by scaling each row of the transition matrix with the exit rate
+            std::vector<typename CtmcType::ValueType>& exitRates = ctmc.getExitRateVector();
+            storm::storage::SparseMatrix<typename CtmcType::ValueType> matrix(std::move(ctmc.getTransitionMatrix()));
+            auto exitRateIt = exitRates.begin();
+            for (uint_fast64_t state = 0; state < matrix.getRowCount(); ++state) {
+                for (auto& entry : matrix.getRow(state)) {
+                    entry.setValue(entry.getValue() / *exitRateIt);
+                }
+                ++exitRateIt;
+            }
+            STORM_LOG_ASSERT(exitRateIt == exitRates.end(), "Unexpected size of rate vector.");
+            
+            // Transform the reward models
+            std::unordered_map<std::string, typename CtmcType::RewardModelType> rewardModels(std::move(ctmc.getRewardModels()));
+            for (auto& rewardModel : rewardModels) {
+                if (rewardModel.second.hasStateRewards()) {
+                    storm::utility::vector::divideVectorsPointwise(rewardModel.second.getStateRewardVector(), exitRates, rewardModel.second.getStateRewardVector());
+                }
+            }
+            
+            if (timeRewardModelName) {
+                // Invert the exit rate vector in place
+                storm::utility::vector::applyPointwise<typename CtmcType::ValueType, typename CtmcType::ValueType>(exitRates, exitRates, [&] (typename CtmcType::ValueType const& r) -> typename CtmcType::ValueType { return storm::utility::one<typename CtmcType::ValueType>() / r; });
+                typename CtmcType::RewardModelType timeRewards(std::move(exitRates));
+                auto insertRes = rewardModels.insert(std::make_pair(*timeRewardModelName, std::move(timeRewards)));
+                STORM_LOG_THROW(insertRes.second, storm::exceptions::InvalidArgumentException, "Could not insert auxiliary reward model " << *timeRewardModelName << " because a model with this name already exists.");
+            }
+            // exitRates might be invalidated at this point!!
+            
+            return std::make_shared<storm::models::sparse::Dtmc<typename CtmcType::ValueType, typename CtmcType::RewardModelType>>(std::move(matrix), std::move(ctmc.getStateLabeling()), std::move(rewardModels));
+        }
+        
+        template<typename CtmcType>
+        bool SparseCtmcToSparseDtmcTransformer<CtmcType>::transformationPreservesProperty(storm::logic::Formula const& formula) {
+            storm::logic::FragmentSpecification fragment = storm::logic::propositional();
+            fragment.setProbabilityOperatorsAllowed(true);
+            fragment.setGloballyFormulasAllowed(true);
+            fragment.setReachabilityProbabilityFormulasAllowed(true);
+            fragment.setNextFormulasAllowed(true);
+            fragment.setUntilFormulasAllowed(true);
+            fragment.setRewardOperatorsAllowed(true);
+            fragment.setReachabilityRewardFormulasAllowed(true);
+            return formula.isInFragment(fragment);
+        }
+        
+        template<typename MaType>
+        std::shared_ptr<storm::models::sparse::Mdp<typename MaType::ValueType, typename MaType::RewardModelType>> SparseMaToSparseMdpTransformer<MaType>::translate(MaType&& ma, boost::optional<std::string> const& timeRewardModelName) {
+            STORM_LOG_THROW(ma.isClosed(), storm::exceptions::InvalidArgumentException, "Transformation of MA to its underlying MDP is only possible for closed MAs");
+            std::vector<typename MaType::ValueType>& exitRates = ma.getExitRates();
+            
+            // Markov automata already store the probability matrix
+            
+            // Transform the reward models
+            std::unordered_map<std::string, typename MaType::RewardModelType> rewardModels(std::move(ma.getRewardModels()));
+            for (auto& rewardModel : rewardModels) {
+                if (rewardModel.second.hasStateRewards()) {
+                    auto& stateRewards = rewardModel.second.getStateRewardVector();
+                    for (auto state : ma.getMarkovianStates()) {
+                        stateRewards[state] /= exitRates[state];
+                    }
+                }
+            }
+            
+            if (timeRewardModelName) {
+                // Invert the exit rate vector. Avoid division by zero at probabilistic states
+                std::vector<typename MaType::ValueType> timeRewardVector(exitRates.size(), storm::utility::zero<typename MaType::ValueType>());
+                for (auto state : ma.getMarkovianStates()) {
+                    timeRewardVector[state] = storm::utility::one<typename MaType::ValueType>() / exitRates[state];
+                }
+                typename MaType::RewardModelType timeRewards(std::move(timeRewardVector));
+                auto insertRes = rewardModels.insert(std::make_pair(*timeRewardModelName, std::move(timeRewards)));
+                STORM_LOG_THROW(insertRes.second, storm::exceptions::InvalidArgumentException, "Could not insert auxiliary reward model " << *timeRewardModelName << " because a model with this name already exists.");
+            }
+            
+            return std::make_shared<storm::models::sparse::Mdp<typename MaType::ValueType, typename MaType::RewardModelType>>(std::move(ma.getTransitionMatrix()), std::move(ma.getStateLabeling()), std::move(rewardModels));
+        }
+        
+        template<typename MaType>
+        bool SparseMaToSparseMdpTransformer<MaType>::transformationPreservesProperty(storm::logic::Formula const& formula) {
+            storm::logic::FragmentSpecification fragment = storm::logic::propositional();
+            fragment.setProbabilityOperatorsAllowed(true);
+            fragment.setGloballyFormulasAllowed(true);
+            fragment.setReachabilityProbabilityFormulasAllowed(true);
+            fragment.setNextFormulasAllowed(true);
+            fragment.setUntilFormulasAllowed(true);
+            fragment.setRewardOperatorsAllowed(true);
+            fragment.setReachabilityRewardFormulasAllowed(true);
+            
+            return formula.isInFragment(fragment);
+        }
+
+        template void transformContinuousToDiscreteModelInPlace<double>(std::shared_ptr<storm::models::sparse::Model<double>>& markovModel, std::shared_ptr<storm::logic::Formula const>& formula);
+        template void transformContinuousToDiscreteModelInPlace<storm::RationalNumber>(std::shared_ptr<storm::models::sparse::Model<storm::RationalNumber>>& markovModel, std::shared_ptr<storm::logic::Formula const>& formula);
+        template void transformContinuousToDiscreteModelInPlace<storm::RationalFunction>(std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>>& markovModel, std::shared_ptr<storm::logic::Formula const>& formula);
+        template class SparseCtmcToSparseDtmcTransformer<storm::models::sparse::Ctmc<double>>;
+        template class SparseCtmcToSparseDtmcTransformer<storm::models::sparse::Ctmc<storm::RationalNumber>>;
+        template class SparseCtmcToSparseDtmcTransformer<storm::models::sparse::Ctmc<storm::RationalFunction>>;
+        template class SparseMaToSparseMdpTransformer<storm::models::sparse::MarkovAutomaton<double>>;
+        template class SparseMaToSparseMdpTransformer<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
+        template class SparseMaToSparseMdpTransformer<storm::models::sparse::MarkovAutomaton<storm::RationalFunction>>;
+    }
+}

src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h

@@ -0,0 +1,43 @@
+#pragma once
+
+#include <boost/optional.hpp>
+
+#include "storm/models/sparse/Dtmc.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/Ctmc.h"
+#include "storm/models/sparse/MarkovAutomaton.h"
+#include "storm/logic/Formula.h"
+
+namespace storm {
+    namespace transformer {
+
+        // Transforms the given continuous model to a discrete time model IN PLACE (i.e., the continuous model will be invalidated).
+        // If such a transformation does not preserve the given formula, the transformation does not take place.
+        // Moreover, the given formula might be changed (e.g. TimeOperatorFormulas become RewardOperatorFormulas).
+        template <typename ValueType>
+        void transformContinuousToDiscreteModelInPlace(std::shared_ptr<storm::models::sparse::Model<ValueType>>& markovModel, std::shared_ptr<storm::logic::Formula const>& formula);
+        
+        template<typename CtmcType>
+        class SparseCtmcToSparseDtmcTransformer {
+        public:
+            // Transforms the given CTMC to its underlying (aka embedded) DTMC.
+            // A reward model for time is added if a corresponding reward model name is given
+            static std::shared_ptr<storm::models::sparse::Dtmc<typename CtmcType::ValueType, typename CtmcType::RewardModelType>> translate(CtmcType&& ctmc, boost::optional<std::string> const& timeRewardModelName = boost::none);
+            
+            // If this method returns true, the given formula is preserced by the transformation
+            static bool transformationPreservesProperty(storm::logic::Formula const& formula);
+            
+        };
+        
+        template<typename MaType>
+        class SparseMaToSparseMdpTransformer {
+        public:
+            // Transforms the given MA to its underlying (aka embedded) MDP.
+            // A reward model for time is added if a corresponding reward model name is given
+            static std::shared_ptr<storm::models::sparse::Mdp<typename MaType::ValueType, typename MaType::RewardModelType>> translate(MaType&& ma, boost::optional<std::string> const& timeRewardModelName = boost::none);
+            
+            // If this method returns true, the given formula is preserved by the transformation
+            static bool transformationPreservesProperty(storm::logic::Formula const& formula);
+        };
+    }
+}

src/storm/transformer/ParameterLifter.cpp

@@ -22,7 +22,7 @@ namespace storm {
             }
             
             // Stores which entries of the original matrix/vector are non-constant. Entries for non-selected rows/columns are omitted
-            storm::storage::BitVector nonConstMatrixEntries(pMatrix.getEntryCount(), false); // note that this vector is too large if there are non-selected rows/columns
+            storm::storage::BitVector nonConstMatrixEntries(pMatrix.getEntryCount(), false); //this vector has to be resized later
             storm::storage::BitVector nonConstVectorEntries(selectedRows.getNumberOfSetBits(), false);
             // Counters for selected entries in the pMatrix and the pVector
             uint_fast64_t pMatrixEntryCount = 0;
@@ -247,7 +247,7 @@ namespace storm {
             // insert the function and the valuation
             auto insertionRes = collectedFunctions.insert(std::pair<FunctionValuation, ConstantType>(FunctionValuation(std::move(simplifiedFunction), std::move(simplifiedValuation)), storm::utility::one<ConstantType>()));
             if(insertionRes.second) {
-                STORM_LOG_THROW(storm::utility::parametric::isMultiLinearPolynomial(insertionRes.first->first.first), storm::exceptions::NotSupportedException, "Parameter lifting for non-multilinear polynomial " << insertionRes.first->first.first << " is not supported");
+                STORM_LOG_WARN_COND(storm::utility::parametric::isMultiLinearPolynomial(insertionRes.first->first.first), "Parameter lifting for non-multilinear polynomial " << insertionRes.first->first.first << " invoked. This might not be sound...");
             }
             //Note that references to elements of an unordered map remain valid after calling unordered_map::insert.
             return insertionRes.first->second;

src/storm/utility/storm.h

@@ -78,6 +78,8 @@
 #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "storm/modelchecker/results/SymbolicQualitativeCheckResult.h"
 
+#include "storm/transformer/ContinuousToDiscreteTimeModelTransformer.h"
+
 // Headers for counterexample generation.
 #include "storm/counterexamples/MILPMinimalLabelSetGenerator.h"
 #include "storm/counterexamples/SMTMinimalCommandSetGenerator.h"
@@ -318,6 +320,14 @@ namespace storm {
     template<>
     inline void performParameterLifting(std::shared_ptr<storm::models::sparse::Model<storm::RationalFunction>> markovModel, std::shared_ptr<storm::logic::Formula const> const& formula) {
         storm::utility::Stopwatch parameterLiftingStopWatch(true);
+        std::shared_ptr<storm::logic::Formula const> consideredFormula = formula;
+        
+        if (markovModel->isOfType(storm::models::ModelType::Ctmc) || markovModel->isOfType(storm::models::ModelType::MarkovAutomaton)) {
+            STORM_PRINT_AND_LOG("Transforming continuous model to discrete model...");
+            storm::transformer::transformContinuousToDiscreteModelInPlace(markovModel, consideredFormula);
+            STORM_PRINT_AND_LOG(" done!" << std::endl);
+        }
+        
         auto modelParameters = storm::models::sparse::getProbabilityParameters(*markovModel);
         auto rewParameters = storm::models::sparse::getRewardParameters(*markovModel);
         modelParameters.insert(rewParameters.begin(), rewParameters.end());
@@ -327,10 +337,10 @@ namespace storm {
         auto parameterSpace = storm::storage::ParameterRegion<storm::RationalFunction>::parseRegion(parameterSpaceAsString, modelParameters);
         auto refinementThreshold = storm::utility::convertNumber<typename storm::storage::ParameterRegion<storm::RationalFunction>::CoefficientType>(storm::settings::getModule<storm::settings::modules::ParametricSettings>().getRefinementThreshold());
         std::vector<std::pair<storm::storage::ParameterRegion<storm::RationalFunction>, storm::modelchecker::parametric::RegionCheckResult>> result;
+                
+        STORM_PRINT_AND_LOG("Performing parameter lifting for property " << *consideredFormula << " with parameter space " << parameterSpace.toString(true) << " and refinement threshold " << storm::utility::convertNumber<double>(refinementThreshold) << " ..." << std::endl);
         
-        STORM_PRINT_AND_LOG("Performing parameter lifting for property " << *formula << " with parameter space " << parameterSpace.toString(true) << " and refinement threshold " << storm::utility::convertNumber<double>(refinementThreshold) << " ..." << std::endl);
-        
-        storm::modelchecker::CheckTask<storm::logic::Formula, storm::RationalFunction> task(*formula, true);
+        storm::modelchecker::CheckTask<storm::logic::Formula, storm::RationalFunction> task(*consideredFormula, true);
         std::string resultVisualization;
         
         if (markovModel->isOfType(storm::models::ModelType::Dtmc)) {

src/storm/utility/vector.h

@@ -391,6 +391,18 @@ namespace storm {
                 applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::multiplies<>());
             }
             
+            /*!
+             * Divides the two given vectors (pointwise) and writes the result to the target vector.
+             *
+             * @param firstOperand The first operand.
+             * @param secondOperand The second operand
+             * @param target The target vector.
+             */
+            template<class InValueType1, class InValueType2, class OutValueType>
+            void divideVectorsPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
+                applyPointwise<InValueType1, InValueType2, OutValueType>(firstOperand, secondOperand, target, std::divides<>());
+            }
+            
             /*!
              * Multiplies each element of the given vector with the given factor and writes the result into the vector.
              *