diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 18157bb94..2d09ecf8f 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -1361,7 +1361,7 @@ namespace storm { } template - void SparseMatrix::scaleRowsInPlace(std::vector const& factors) { + void SparseMatrix::scaleRowsInPlace(std::vector const& factors) { STORM_LOG_ASSERT(factors.size() == this->getRowCount(), "Can not scale rows: Number of rows and number of scaling factors do not match."); uint_fast64_t row = 0; for (auto const& factor : factors) { @@ -1371,7 +1371,26 @@ namespace storm { ++row; } } - + + template + void SparseMatrix::divideRowsInPlace(std::vector const& divisors) { + STORM_LOG_ASSERT(divisors.size() == this->getRowCount(), "Can not divide rows: Number of rows and number of divisors do not match."); + uint_fast64_t row = 0; + for (auto const& divisor : divisors) { + STORM_LOG_ASSERT(!storm::utility::isZero(divisor), "Can not divide row " << row << " by 0."); + for (auto& entry : getRow(row)) { + entry.setValue(entry.getValue() / divisor); + } + ++row; + } + } + +#ifdef STORM_HAVE_CARL + template<> + void SparseMatrix::divideRowsInPlace(std::vector const&) { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); + } +#endif template typename SparseMatrix::const_rows SparseMatrix::getRows(index_type startRow, index_type endRow) const { diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index 01e4079a0..f3b9f280f 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -796,9 +796,18 @@ s * @param insertDiagonalEntries If set to true, the resulting matri /*! * Scales each row of the matrix, i.e., multiplies each element in row i with factors[i] + * + * @param factors The factors with which each row is scaled. */ void scaleRowsInPlace(std::vector const& factors); + /*! + * Divides each row of the matrix, i.e., divides each element in row i with divisors[i] + * + * @param divisors The divisors with which each row is divided. + */ + void divideRowsInPlace(std::vector const& divisors); + /*! * Performs one step of the successive over-relaxation technique. * diff --git a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp index 4f1058c5c..f45b7076c 100644 --- a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp +++ b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.cpp @@ -14,8 +14,8 @@ namespace storm { namespace transformer { - template - void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula) { + template + std::shared_ptr> transformContinuousToDiscreteModel(std::shared_ptr> markovModel, std::shared_ptr& formula) { boost::optional timeRewardModelName; if (formula->isTimeOperatorFormula()) { auto const& timeOpFormula = formula->asTimeOperatorFormula(); @@ -29,13 +29,42 @@ namespace storm { } if (markovModel->isOfType(storm::models::ModelType::Ctmc)) { - SparseCtmcToSparseDtmcTransformer> transformer; + SparseCtmcToSparseDtmcTransformer transformer; + if (transformer.transformationPreservesProperty(*formula)) { + STORM_LOG_INFO("Transforming Ctmc to embedded Dtmc..."); + markovModel = transformer.translate(*markovModel->template as>(), timeRewardModelName); + } + } else if (markovModel->isOfType(storm::models::ModelType::MarkovAutomaton)) { + SparseMaToSparseMdpTransformer transformer; + if (transformer.transformationPreservesProperty(*formula)) { + STORM_LOG_INFO("Transforming Markov automaton to embedded Mdp..."); + markovModel = transformer.translate(*markovModel->template as>(), timeRewardModelName); + } + } + } + + template + void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula) { + boost::optional timeRewardModelName; + if (formula->isTimeOperatorFormula()) { + auto const& timeOpFormula = formula->asTimeOperatorFormula(); + if (timeOpFormula.getSubformula().isReachabilityTimeFormula()) { + auto reachabilityRewardFormula = std::make_shared(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(reachabilityRewardFormula, timeRewardModelName, timeOpFormula.getOperatorInformation()); + } + } + + if (markovModel->isOfType(storm::models::ModelType::Ctmc)) { + SparseCtmcToSparseDtmcTransformer transformer; if (transformer.transformationPreservesProperty(*formula)) { STORM_LOG_INFO("Transforming Ctmc to embedded Dtmc..."); markovModel = transformer.translate(std::move(*markovModel->template as>()), timeRewardModelName); } } else if (markovModel->isOfType(storm::models::ModelType::MarkovAutomaton)) { - SparseMaToSparseMdpTransformer> transformer; + SparseMaToSparseMdpTransformer transformer; if (transformer.transformationPreservesProperty(*formula)) { STORM_LOG_INFO("Transforming Markov automaton to embedded Mdp..."); markovModel = transformer.translate(std::move(*markovModel->template as>()), timeRewardModelName); @@ -43,23 +72,54 @@ namespace storm { } } - template - std::shared_ptr> SparseCtmcToSparseDtmcTransformer::translate(CtmcType&& ctmc, boost::optional const& timeRewardModelName) { - // Turn the rates into probabilities by scaling each row of the transition matrix with the exit rate - std::vector& exitRates = ctmc.getExitRateVector(); - storm::storage::SparseMatrix 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); + template + std::shared_ptr> SparseCtmcToSparseDtmcTransformer::translate(storm::models::sparse::Ctmc const& ctmc, boost::optional const& timeRewardModelName) { + // Init the dtmc components + storm::storage::sparse::ModelComponents dtmcComponents(ctmc.getTransitionMatrix(), ctmc.getStateLabeling(), ctmc.getRewardModels()); + dtmcComponents.choiceLabeling = ctmc.getOptionalChoiceLabeling(); + dtmcComponents.stateValuations = ctmc.getOptionalStateValuations(); + dtmcComponents.choiceOrigins = ctmc.getOptionalChoiceOrigins(); + + // Turn the rates into probabilities by dividing each row of the transition matrix with the exit rate + std::vector const& exitRates = ctmc.getExitRateVector(); + dtmcComponents.transitionMatrix.divideRowsInPlace(exitRates); + + // Transform the reward models + for (auto& rewardModel : dtmcComponents.rewardModels) { + if (rewardModel.second.hasStateRewards()) { + storm::utility::vector::divideVectorsPointwise(rewardModel.second.getStateRewardVector(), exitRates, rewardModel.second.getStateRewardVector()); } - ++exitRateIt; } - STORM_LOG_ASSERT(exitRateIt == exitRates.end(), "Unexpected size of rate vector."); + + if (timeRewardModelName) { + // Invert the exit rate vector in place + std::vector timeRewardVector; + timeRewardVector.reserve(exitRates.size()); + for (auto const& r : exitRates) { + timeRewardVector.push_back(storm::utility::one() / r); + } + RewardModelType timeRewards(std::move(timeRewardVector)); + auto insertRes = dtmcComponents.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>(std::move(dtmcComponents)); + } + + template + std::shared_ptr> SparseCtmcToSparseDtmcTransformer::translate(storm::models::sparse::Ctmc&& ctmc, boost::optional const& timeRewardModelName) { + // Init the dtmc components + storm::storage::sparse::ModelComponents dtmcComponents(std::move(ctmc.getTransitionMatrix()), std::move(ctmc.getStateLabeling()), std::move(ctmc.getRewardModels())); + dtmcComponents.choiceLabeling = std::move(ctmc.getOptionalChoiceLabeling()); + dtmcComponents.stateValuations = std::move(ctmc.getOptionalStateValuations()); + dtmcComponents.choiceOrigins = std::move(ctmc.getOptionalChoiceOrigins()); + + // Turn the rates into probabilities by dividing each row of the transition matrix with the exit rate + std::vector& exitRates = ctmc.getExitRateVector(); + dtmcComponents.transitionMatrix.divideRowsInPlace(exitRates); // Transform the reward models - std::unordered_map rewardModels(std::move(ctmc.getRewardModels())); - for (auto& rewardModel : rewardModels) { + for (auto& rewardModel : dtmcComponents.rewardModels) { if (rewardModel.second.hasStateRewards()) { storm::utility::vector::divideVectorsPointwise(rewardModel.second.getStateRewardVector(), exitRates, rewardModel.second.getStateRewardVector()); } @@ -67,18 +127,18 @@ namespace storm { if (timeRewardModelName) { // Invert the exit rate vector in place - storm::utility::vector::applyPointwise(exitRates, exitRates, [&] (typename CtmcType::ValueType const& r) -> typename CtmcType::ValueType { return storm::utility::one() / r; }); - typename CtmcType::RewardModelType timeRewards(std::move(exitRates)); - auto insertRes = rewardModels.insert(std::make_pair(*timeRewardModelName, std::move(timeRewards))); + storm::utility::vector::applyPointwise(exitRates, exitRates, [&] (ValueType const& r) -> ValueType { return storm::utility::one() / r; }); + RewardModelType timeRewards(std::move(exitRates)); + auto insertRes = dtmcComponents.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!! + // Note: exitRates might be invalidated at this point. - return std::make_shared>(std::move(matrix), std::move(ctmc.getStateLabeling()), std::move(rewardModels)); + return std::make_shared>(std::move(dtmcComponents)); } - - template - bool SparseCtmcToSparseDtmcTransformer::transformationPreservesProperty(storm::logic::Formula const& formula) { + + template + bool SparseCtmcToSparseDtmcTransformer::transformationPreservesProperty(storm::logic::Formula const& formula) { storm::logic::FragmentSpecification fragment = storm::logic::propositional(); fragment.setProbabilityOperatorsAllowed(true); fragment.setGloballyFormulasAllowed(true); @@ -90,16 +150,22 @@ namespace storm { return formula.isInFragment(fragment); } - template - std::shared_ptr> SparseMaToSparseMdpTransformer::translate(MaType&& ma, boost::optional const& timeRewardModelName) { + + template + std::shared_ptr> SparseMaToSparseMdpTransformer::translate(storm::models::sparse::MarkovAutomaton const& ma, boost::optional 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& exitRates = ma.getExitRates(); + + // Init the mdp components + storm::storage::sparse::ModelComponents mdpComponents(ma.getTransitionMatrix(), ma.getStateLabeling(), ma.getRewardModels()); + mdpComponents.choiceLabeling = ma.getOptionalChoiceLabeling(); + mdpComponents.stateValuations = ma.getOptionalStateValuations(); + mdpComponents.choiceOrigins = ma.getOptionalChoiceOrigins(); // Markov automata already store the probability matrix // Transform the reward models - std::unordered_map rewardModels(std::move(ma.getRewardModels())); - for (auto& rewardModel : rewardModels) { + std::vector const& exitRates = ma.getExitRates(); + for (auto& rewardModel : mdpComponents.rewardModels) { if (rewardModel.second.hasStateRewards()) { auto& stateRewards = rewardModel.second.getStateRewardVector(); for (auto state : ma.getMarkovianStates()) { @@ -110,20 +176,57 @@ namespace storm { if (timeRewardModelName) { // Invert the exit rate vector. Avoid division by zero at probabilistic states - std::vector timeRewardVector(exitRates.size(), storm::utility::zero()); + std::vector timeRewardVector(exitRates.size(), storm::utility::zero()); for (auto state : ma.getMarkovianStates()) { - timeRewardVector[state] = storm::utility::one() / exitRates[state]; + timeRewardVector[state] = storm::utility::one() / exitRates[state]; } - typename MaType::RewardModelType timeRewards(std::move(timeRewardVector)); - auto insertRes = rewardModels.insert(std::make_pair(*timeRewardModelName, std::move(timeRewards))); + RewardModelType timeRewards(std::move(timeRewardVector)); + auto insertRes = mdpComponents.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>(std::move(ma.getTransitionMatrix()), std::move(ma.getStateLabeling()), std::move(rewardModels)); + return std::make_shared>(std::move(mdpComponents)); + } + + template + std::shared_ptr> SparseMaToSparseMdpTransformer::translate(storm::models::sparse::MarkovAutomaton&& ma, boost::optional 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& exitRates = ma.getExitRates(); + + // Init the mdp components + storm::storage::sparse::ModelComponents mdpComponents(std::move(ma.getTransitionMatrix()), std::move(ma.getStateLabeling()), std::move(ma.getRewardModels())); + mdpComponents.choiceLabeling = std::move(ma.getOptionalChoiceLabeling()); + mdpComponents.stateValuations = std::move(ma.getOptionalStateValuations()); + mdpComponents.choiceOrigins = std::move(ma.getOptionalChoiceOrigins()); + + // Markov automata already store the probability matrix + + // Transform the reward models + for (auto& rewardModel : mdpComponents.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 timeRewardVector(exitRates.size(), storm::utility::zero()); + for (auto state : ma.getMarkovianStates()) { + timeRewardVector[state] = storm::utility::one() / exitRates[state]; + } + RewardModelType timeRewards(std::move(timeRewardVector)); + auto insertRes = mdpComponents.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>(std::move(mdpComponents)); } - template - bool SparseMaToSparseMdpTransformer::transformationPreservesProperty(storm::logic::Formula const& formula) { + template + bool SparseMaToSparseMdpTransformer::transformationPreservesProperty(storm::logic::Formula const& formula) { storm::logic::FragmentSpecification fragment = storm::logic::propositional(); fragment.setProbabilityOperatorsAllowed(true); fragment.setGloballyFormulasAllowed(true); @@ -135,15 +238,18 @@ namespace storm { return formula.isInFragment(fragment); } - + + template std::shared_ptr> transformContinuousToDiscreteModel(std::shared_ptr> markovModel, std::shared_ptr& formula); + template std::shared_ptr> transformContinuousToDiscreteModel(std::shared_ptr> markovModel, std::shared_ptr& formula); + template std::shared_ptr> transformContinuousToDiscreteModel(std::shared_ptr> markovModel, std::shared_ptr& formula); template void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula); template void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula); template void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula); - template class SparseCtmcToSparseDtmcTransformer>; - template class SparseCtmcToSparseDtmcTransformer>; - template class SparseCtmcToSparseDtmcTransformer>; - template class SparseMaToSparseMdpTransformer>; - template class SparseMaToSparseMdpTransformer>; - template class SparseMaToSparseMdpTransformer>; + template class SparseCtmcToSparseDtmcTransformer; + template class SparseCtmcToSparseDtmcTransformer; + template class SparseCtmcToSparseDtmcTransformer; + template class SparseMaToSparseMdpTransformer; + template class SparseMaToSparseMdpTransformer; + template class SparseMaToSparseMdpTransformer; } } diff --git a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h index e3985924c..cea1bf367 100644 --- a/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h +++ b/src/storm/transformer/ContinuousToDiscreteTimeModelTransformer.h @@ -6,35 +6,44 @@ #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/Ctmc.h" #include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/models/sparse/StandardRewardModel.h" #include "storm/logic/Formula.h" namespace storm { namespace transformer { + + // Transforms the given continuous model to a discrete time model. + // If such a transformation does not preserve the given formula, the transformation does not take place and the original model is returned + // Moreover, the given formula might be changed (e.g. TimeOperatorFormulas become RewardOperatorFormulas). + template > + std::shared_ptr> transformContinuousToDiscreteModel(std::shared_ptr> markovModel, std::shared_ptr& formula); + // 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 - void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula); + template > + void transformContinuousToDiscreteModelInPlace(std::shared_ptr>& markovModel, std::shared_ptr& formula); - template + template > 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> translate(CtmcType&& ctmc, boost::optional const& timeRewardModelName = boost::none); + static std::shared_ptr> translate(storm::models::sparse::Ctmc const& ctmc, boost::optional const& timeRewardModelName = boost::none); + static std::shared_ptr> translate(storm::models::sparse::Ctmc&& ctmc, boost::optional 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 + template > 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> translate(MaType&& ma, boost::optional const& timeRewardModelName = boost::none); + static std::shared_ptr> translate(storm::models::sparse::MarkovAutomaton const& ma, boost::optional const& timeRewardModelName = boost::none); + static std::shared_ptr> translate(storm::models::sparse::MarkovAutomaton&& ma, boost::optional 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);