diff --git a/src/storm/modelchecker/parametric/ParameterLifting.cpp b/src/storm/modelchecker/parametric/ParameterLifting.cpp index 4db278bba..59704a120 100644 --- a/src/storm/modelchecker/parametric/ParameterLifting.cpp +++ b/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(); diff --git a/src/storm/models/sparse/Ctmc.cpp b/src/storm/models/sparse/Ctmc.cpp index 884eb359d..c2913a254 100644 --- a/src/storm/models/sparse/Ctmc.cpp +++ b/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()); diff --git a/src/storm/models/sparse/Ctmc.h b/src/storm/models/sparse/Ctmc.h index 417a899a4..3ebd5249a 100644 --- a/src/storm/models/sparse/Ctmc.h +++ b/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: /*! diff --git a/src/storm/storage/ParameterRegion.cpp b/src/storm/storage/ParameterRegion.cpp index 3f71db902..9e5287d05 100644 --- a/src/storm/storage/ParameterRegion.cpp +++ b/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)); diff --git a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp new file mode 100644 index 000000000..4f1058c5c --- /dev/null +++ b/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>>; + } +} diff --git a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h new file mode 100644 index 000000000..e3985924c --- /dev/null +++ b/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); + }; + } +} diff --git a/src/storm/transformer/ParameterLifter.cpp b/src/storm/transformer/ParameterLifter.cpp index 9fe1d333d..ffd2d1bd5 100644 --- a/src/storm/transformer/ParameterLifter.cpp +++ b/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; diff --git a/src/storm/utility/storm.h b/src/storm/utility/storm.h index eb6dcf52d..7749495e2 100644 --- a/src/storm/utility/storm.h +++ b/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)) { diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 578a942a8..77932dab6 100644 --- a/src/storm/utility/vector.h +++ b/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. *