From a1ad5377e83e14733b1455e3a7409f32d707b980 Mon Sep 17 00:00:00 2001 From: TimQu Date: Mon, 6 Mar 2017 13:54:21 +0100 Subject: [PATCH] implemented parameter lifting model checker for DTMCs (as a replacement of the old ApproximationModel). --- src/storm/modelchecker/CheckTask.h | 15 ++ .../parametric/ParameterRegion.cpp | 251 ++++++++++++++++++ .../modelchecker/parametric/ParameterRegion.h | 114 ++++++++ ...SparseDtmcParameterLiftingModelChecker.cpp | 201 ++++++++++++++ .../SparseDtmcParameterLiftingModelChecker.h | 47 ++++ .../SparseInstantiationModelChecker.cpp | 3 +- .../SparseParameterLiftingModelChecker.cpp | 92 +++++++ .../SparseParameterLiftingModelChecker.h | 55 ++++ .../region/SparseMdpRegionModelChecker.cpp | 2 +- src/storm/transformer/ParameterLifter.cpp | 248 +++++++++++++++++ src/storm/transformer/ParameterLifter.h | 133 ++++++++++ .../SparseParametricModelSimplifier.cpp | 3 - src/storm/utility/ModelInstantiator.cpp | 7 - src/storm/utility/ModelInstantiator.h | 7 - src/storm/utility/constants.cpp | 5 + src/storm/utility/parametric.cpp | 12 +- src/storm/utility/parametric.h | 8 +- 17 files changed, 1175 insertions(+), 28 deletions(-) create mode 100644 src/storm/modelchecker/parametric/ParameterRegion.cpp create mode 100644 src/storm/modelchecker/parametric/ParameterRegion.h create mode 100644 src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.cpp create mode 100644 src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.h create mode 100644 src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.cpp create mode 100644 src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h create mode 100644 src/storm/transformer/ParameterLifter.cpp create mode 100644 src/storm/transformer/ParameterLifter.h diff --git a/src/storm/modelchecker/CheckTask.h b/src/storm/modelchecker/CheckTask.h index 3caba461c..7c1249e11 100644 --- a/src/storm/modelchecker/CheckTask.h +++ b/src/storm/modelchecker/CheckTask.h @@ -3,6 +3,7 @@ #include #include +#include #include "storm/logic/Formulas.h" #include "storm/utility/constants.h" @@ -84,6 +85,20 @@ namespace storm { CheckTask substituteFormula(NewFormulaType const& newFormula) const { return CheckTask(newFormula, this->optimizationDirection, this->rewardModel, this->onlyInitialStatesRelevant, this->bound, this->qualitative, this->produceSchedulers, this->hint); } + + /*! + * Copies the check task from the source while replacing the considered ValueType the new one. In particular, this + * changes the formula type of the check task object. + */ + template + CheckTask convertValueType() const { + boost::optional> newBound; + if(this->bound.is_initialized()) { + STORM_LOG_THROW(storm::utility::isConstant(this->getBoundThreshold()), storm::exceptions::InvalidTypeException, "Tried to convert a non-constant threshold "); + newBound = storm::logic::Bound(this->getBoundComparisonType(), storm::utility::convertNumber(this->getBoundThreshold())); + } + return CheckTask(this->formula, this->optimizationDirection, this->rewardModel, this->onlyInitialStatesRelevant, newBound, this->qualitative, this->produceSchedulers, this->hint); + } /*! * Retrieves the formula from this task. diff --git a/src/storm/modelchecker/parametric/ParameterRegion.cpp b/src/storm/modelchecker/parametric/ParameterRegion.cpp new file mode 100644 index 000000000..aab7e9d06 --- /dev/null +++ b/src/storm/modelchecker/parametric/ParameterRegion.cpp @@ -0,0 +1,251 @@ +#include "storm/modelchecker/parametric/ParameterRegion.h" + +#include "storm/utility/region.h" +#include "storm/utility/macros.h" +#include "storm/parser/MappedFile.h" +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/RegionSettings.h" +#include "storm/exceptions/InvalidSettingsException.h" +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/utility/constants.h" +#include "storm/utility/file.h" + +namespace storm { + namespace modelchecker { + namespace parametric { + + template + ParameterRegion::ParameterRegion(Valuation const& lowerBoundaries, Valuation const& upperBoundaries) : lowerBoundaries(lowerBoundaries), upperBoundaries(upperBoundaries) { + init(); + } + + template + ParameterRegion::ParameterRegion(Valuation&& lowerBoundaries, Valuation&& upperBoundaries) : lowerBoundaries(lowerBoundaries), upperBoundaries(upperBoundaries) { + init(); + } + + template + void ParameterRegion::init() { + //check whether both mappings map the same variables, check that lower boundary <= upper boundary, and pre-compute the set of variables + for (auto const& variableWithLowerBoundary : this->lowerBoundaries) { + auto variableWithUpperBoundary = this->upperBoundaries.find(variableWithLowerBoundary.first); + STORM_LOG_THROW((variableWithUpperBoundary != upperBoundaries.end()), storm::exceptions::InvalidArgumentException, "Could not create region. No upper boundary specified for Variable " << variableWithLowerBoundary.first); + STORM_LOG_THROW((variableWithLowerBoundary.second<=variableWithUpperBoundary->second), storm::exceptions::InvalidArgumentException, "Could not create region. The lower boundary for " << variableWithLowerBoundary.first << " is larger then the upper boundary"); + this->variables.insert(variableWithLowerBoundary.first); + } + for (auto const& variableWithBoundary : this->upperBoundaries) { + STORM_LOG_THROW((this->variables.find(variableWithBoundary.first) != this->variables.end()), storm::exceptions::InvalidArgumentException, "Could not create region. No lower boundary specified for Variable " << variableWithBoundary.first); + } + } + + template + std::set::VariableType> const& ParameterRegion::getVariables() const { + return this->variables; + } + + template + typename ParameterRegion::CoefficientType const& ParameterRegion::getLowerBoundary(VariableType const& variable) const { + auto const& result = lowerBoundaries.find(variable); + STORM_LOG_THROW(result != lowerBoundaries.end(), storm::exceptions::InvalidArgumentException, "Tried to find a lower boundary for variable " << variable << " which is not specified by this region"); + return (*result).second; + } + + template + typename ParameterRegion::CoefficientType const& ParameterRegion::getUpperBoundary(VariableType const& variable) const { + auto const& result = upperBoundaries.find(variable); + STORM_LOG_THROW(result != upperBoundaries.end(), storm::exceptions::InvalidArgumentException, "Tried to find an upper boundary for variable " << variable << " which is not specified by this region"); + return (*result).second; + } + + template + typename ParameterRegion::Valuation const& ParameterRegion::getUpperBoundaries() const { + return upperBoundaries; + } + + template + typename ParameterRegion::Valuation const& ParameterRegion::getLowerBoundaries() const { + return lowerBoundaries; + } + + template + std::vector::Valuation> ParameterRegion::getVerticesOfRegion(std::set const& consideredVariables) const { + std::size_t const numOfVariables = consideredVariables.size(); + std::size_t const numOfVertices = std::pow(2, numOfVariables); + std::vector resultingVector(numOfVertices); + + for (uint_fast64_t vertexId = 0; vertexId < numOfVertices; ++vertexId) { + //interprete vertexId as a bit sequence + //the consideredVariables.size() least significant bits of vertex will always represent the next vertex + //(00...0 = lower boundaries for all variables, 11...1 = upper boundaries for all variables) + uint_fast64_t variableIndex = 0; + for (auto const& variable : consideredVariables) { + if ((vertexId >> variableIndex) % 2 == 0) { + resultingVector[vertexId].insert(std::pair(variable, getLowerBoundary(variable))); + } else { + resultingVector[vertexId].insert(std::pair(variable, getUpperBoundary(variable))); + } + ++variableIndex; + } + } + return resultingVector; + } + + template + typename ParameterRegion::Valuation ParameterRegion::getSomePoint() const { + return this->getLowerBoundaries(); + } + + template + typename ParameterRegion::Valuation ParameterRegion::getCenterPoint() const { + Valuation result; + for (auto const& variable : this->variables) { + result.insert(typename Valuation::value_type(variable, (this->getLowerBoundary(variable) + this->getUpperBoundary(variable))/2)); + } + return result; + } + + template + typename ParameterRegion::CoefficientType ParameterRegion::area() const { + CoefficientType result = storm::utility::one(); + for( auto const& variable : this->variables){ + result *= (this->getUpperBoundary(variable) - this->getLowerBoundary(variable)); + } + return result; + } + + template + void ParameterRegion::split(Valuation const& splittingPoint, std::vector >& regionVector) const{ + //Check if splittingPoint is valid. + STORM_LOG_THROW(splittingPoint.size() == this->variables.size(), storm::exceptions::InvalidArgumentException, "Tried to split a region w.r.t. a point, but the point considers a different number of variables."); + for(auto const& variable : this->variables){ + auto splittingPointEntry=splittingPoint.find(variable); + STORM_LOG_THROW(splittingPointEntry != splittingPoint.end(), storm::exceptions::InvalidArgumentException, "Tried to split a region but a variable of this region is not defined by the splitting point."); + STORM_LOG_THROW(this->getLowerBoundary(variable) <=splittingPointEntry->second, storm::exceptions::InvalidArgumentException, "Tried to split a region but the splitting point is not contained in the region."); + STORM_LOG_THROW(this->getUpperBoundary(variable) >=splittingPointEntry->second, storm::exceptions::InvalidArgumentException, "Tried to split a region but the splitting point is not contained in the region."); + } + + //Now compute the subregions. + std::vector vertices(this->getVerticesOfRegion(this->variables)); + for(auto const& vertex : vertices){ + //The resulting subregion is the smallest region containing vertex and splittingPoint. + Valuation subLower, subUpper; + for(auto variableBound : this->lowerBoundaries){ + VariableType variable = variableBound.first; + auto vertexEntry=vertex.find(variable); + auto splittingPointEntry=splittingPoint.find(variable); + subLower.insert(typename Valuation::value_type(variable, std::min(vertexEntry->second, splittingPointEntry->second))); + subUpper.insert(typename Valuation::value_type(variable, std::max(vertexEntry->second, splittingPointEntry->second))); + } + ParameterRegion subRegion(std::move(subLower), std::move(subUpper)); + if(!storm::utility::isZero(subRegion.area())){ + regionVector.push_back(std::move(subRegion)); + } + } + } + + template + std::string ParameterRegion::toString(bool boundariesAsDouble) const { + std::stringstream regionstringstream; + if(boundariesAsDouble) { + for (auto var : this->getVariables()) { + regionstringstream << storm::utility::convertNumber(this->getLowerBoundary(var)); + regionstringstream << "<="; + regionstringstream << var; + regionstringstream << "<="; + regionstringstream << storm::utility::convertNumber(this->getUpperBoundary(var)); + regionstringstream << ","; + } + } else { + for (auto var : this->getVariables()) { + regionstringstream << this->getLowerBoundary(var); + regionstringstream << "<="; + regionstringstream << var; + regionstringstream << "<="; + regionstringstream << this->getUpperBoundary(var); + regionstringstream << ","; + } + } + std::string regionstring = regionstringstream.str(); + //the last comma should actually be a semicolon + regionstring = regionstring.substr(0, regionstring.length() - 1) + ";"; + return regionstring; + } + + + template + void ParameterRegion::parseParameterBoundaries( + Valuation& lowerBoundaries, + Valuation& upperBoundaries, + std::string const& parameterBoundariesString){ + + std::string::size_type positionOfFirstRelation = parameterBoundariesString.find("<="); + STORM_LOG_THROW(positionOfFirstRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a '<=' after the first number"); + std::string::size_type positionOfSecondRelation = parameterBoundariesString.find("<=", positionOfFirstRelation+2); + STORM_LOG_THROW(positionOfSecondRelation!=std::string::npos, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a '<=' after the parameter"); + + std::string parameter=parameterBoundariesString.substr(positionOfFirstRelation+2,positionOfSecondRelation-(positionOfFirstRelation+2)); + //removes all whitespaces from the parameter string: + parameter.erase(std::remove_if(parameter.begin(), parameter.end(), ::isspace), parameter.end()); + STORM_LOG_THROW(parameter.length()>0, storm::exceptions::InvalidArgumentException, "When parsing the region" << parameterBoundariesString << " I could not find a parameter"); + + VariableType var = storm::utility::region::getVariableFromString(parameter); + CoefficientType lb = storm::utility::convertNumber(parameterBoundariesString.substr(0,positionOfFirstRelation)); + CoefficientType ub = storm::utility::convertNumber(parameterBoundariesString.substr(positionOfSecondRelation+2)); + lowerBoundaries.emplace(std::make_pair(var, lb)); + upperBoundaries.emplace(std::make_pair(var, ub)); + } + + template + ParameterRegion ParameterRegion::parseRegion( + std::string const& regionString){ + Valuation lowerBoundaries; + Valuation upperBoundaries; + std::vector parameterBoundaries; + boost::split(parameterBoundaries, regionString, boost::is_any_of(",")); + for(auto const& parameterBoundary : parameterBoundaries){ + if(!std::all_of(parameterBoundary.begin(),parameterBoundary.end(), ::isspace)){ //skip this string if it only consists of space + parseParameterBoundaries(lowerBoundaries, upperBoundaries, parameterBoundary); + } + } + return ParameterRegion(std::move(lowerBoundaries), std::move(upperBoundaries)); + } + + template + std::vector> ParameterRegion::parseMultipleRegions( + std::string const& regionsString){ + std::vector result; + std::vector regionsStrVec; + boost::split(regionsStrVec, regionsString, boost::is_any_of(";")); + for(auto const& regionStr : regionsStrVec){ + if(!std::all_of(regionStr.begin(),regionStr.end(), ::isspace)){ //skip this string if it only consists of space + result.emplace_back(parseRegion(regionStr)); + } + } + return result; + } + + template + std::vector> ParameterRegion::getRegionsFromSettings(){ + STORM_LOG_THROW(storm::settings::getModule().isRegionsSet() ||storm::settings::getModule().isRegionFileSet(), storm::exceptions::InvalidSettingsException, "Tried to obtain regions from the settings but no regions are specified."); + STORM_LOG_THROW(!(storm::settings::getModule().isRegionsSet() && storm::settings::getModule().isRegionFileSet()), storm::exceptions::InvalidSettingsException, "Regions are specified via file AND cmd line. Only one option is allowed."); + + std::string regionsString; + if(storm::settings::getModule().isRegionsSet()){ + regionsString = storm::settings::getModule().getRegionsFromCmdLine(); + } + else{ + //if we reach this point we can assume that the region is given as a file. + STORM_LOG_THROW(storm::utility::fileExistsAndIsReadable(storm::settings::getModule().getRegionFilePath()), storm::exceptions::InvalidSettingsException, "The path to the file in which the regions are specified is not valid."); + storm::parser::MappedFile mf(storm::settings::getModule().getRegionFilePath().c_str()); + regionsString = std::string(mf.getData(),mf.getDataSize()); + } + return parseMultipleRegions(regionsString); + } + +#ifdef STORM_HAVE_CARL + template class ParameterRegion; +#endif + } + } +} + diff --git a/src/storm/modelchecker/parametric/ParameterRegion.h b/src/storm/modelchecker/parametric/ParameterRegion.h new file mode 100644 index 000000000..2c07abbdb --- /dev/null +++ b/src/storm/modelchecker/parametric/ParameterRegion.h @@ -0,0 +1,114 @@ +#pragma once + +#include + +#include "storm/utility/parametric.h" + +namespace storm { + namespace modelchecker{ + namespace parametric { + template + class ParameterRegion{ + public: + typedef typename storm::utility::parametric::VariableType::type VariableType; + typedef typename storm::utility::parametric::CoefficientType::type CoefficientType; + typedef typename storm::utility::parametric::Valuation Valuation; + + ParameterRegion(Valuation const& lowerBoundaries, Valuation const& upperBoundaries); + ParameterRegion(Valuation&& lowerBoundaries, Valuation&& upperBoundaries); + ParameterRegion(ParameterRegion const& other) = default; + ParameterRegion(ParameterRegion&& other) = default; + + virtual ~ParameterRegion() = default; + + std::set const& getVariables() const; + CoefficientType const& getLowerBoundary(VariableType const& variable) const; + CoefficientType const& getUpperBoundary(VariableType const& variable) const; + Valuation const& getLowerBoundaries() const; + Valuation const& getUpperBoundaries() const; + + /*! + * Returns a vector of all possible combinations of lower and upper bounds of the given variables. + * The first entry of the returned vector will map every variable to its lower bound + * The second entry will map every variable to its lower bound, except the first one (i.e. *getVariables.begin()) + * ... + * The last entry will map every variable to its upper bound + * + * If the given set of variables is empty, the returned vector will contain an empty map + */ + std::vector getVerticesOfRegion(std::set const& consideredVariables) const; + + /*! + * Returns some point that lies within this region + */ + Valuation getSomePoint() const; + + /*! + * Returns the center point of this region + */ + Valuation getCenterPoint() const; + + /*! + * Returns the area of this region + */ + CoefficientType area() const; + + /*! + * Splits the region at the given point and inserts the resulting subregions at the end of the given vector. + * It is assumed that the point lies within this region. + * Subregions with area()==0 are not inserted in the vector. + */ + void split(Valuation const& splittingPoint, std::vector>& regionVector) const; + + //returns the region as string in the format 0.3<=p<=0.4,0.2<=q<=0.5; + std::string toString(bool boundariesAsDouble = false) const; + + /* + * Can be used to parse a single parameter with its boundaries from a string of the form "0.3<=p<=0.5". + * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType. + * The results will be inserted in the given maps + * + */ + static void parseParameterBoundaries( + Valuation& lowerBoundaries, + Valuation& upperBoundaries, + std::string const& parameterBoundariesString + ); + + /* + * Can be used to parse a single region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7". + * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType. + * + */ + static ParameterRegion parseRegion( + std::string const& regionString + ); + + /* + * Can be used to parse a vector of region from a string of the form "0.3<=p<=0.5,0.4<=q<=0.7;0.1<=p<=0.3,0.2<=q<=0.4". + * The numbers are parsed as doubles and then converted to SparseDtmcRegionModelChecker::CoefficientType. + * + */ + static std::vector parseMultipleRegions( + std::string const& regionsString + ); + + + /* + * Retrieves the regions that are specified in the settings. + */ + static std::vector getRegionsFromSettings(); + + private: + + void init(); + + Valuation lowerBoundaries; + Valuation upperBoundaries; + std::set variables; + }; + } + } +} + + diff --git a/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.cpp b/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.cpp new file mode 100644 index 000000000..fd613dcf6 --- /dev/null +++ b/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.cpp @@ -0,0 +1,201 @@ +#include "SparseDtmcParameterLiftingModelChecker.h" + +#include "storm/adapters/CarlAdapter.h" +#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/models/sparse/Dtmc.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/utility/vector.h" +#include "storm/utility/graph.h" + +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/exceptions/NotSupportedException.h" + +namespace storm { + namespace modelchecker { + namespace parametric { + + template + SparseDtmcParameterLiftingModelChecker::SparseDtmcParameterLiftingModelChecker(SparseModelType const& parametricModel) : SparseParameterLiftingModelChecker(parametricModel), solverFactory(std::make_unique>()) { + solverFactory->setTrackScheduler(true); + } + + template + SparseDtmcParameterLiftingModelChecker::SparseDtmcParameterLiftingModelChecker(SparseModelType const& parametricModel, std::unique_ptr>&& solverFactory) : SparseParameterLiftingModelChecker(this->parametricModel), solverFactory(std::move(solverFactory)) { + solverFactory->setTrackScheduler(true); + } + + template + void SparseDtmcParameterLiftingModelChecker::specifyBoundedUntilFormula(CheckTask const& checkTask) { + + // get the step bound + STORM_LOG_THROW(!checkTask.getFormula().hasLowerBound(), storm::exceptions::NotSupportedException, "Lower step bounds are not supported."); + STORM_LOG_THROW(checkTask.getFormula().hasUpperBound(), storm::exceptions::NotSupportedException, "Expected a bounded until formula with an upper bound."); + STORM_LOG_THROW(checkTask.getFormula().isStepBounded(), storm::exceptions::NotSupportedException, "Expected a bounded until formula with step bounds."); + stepBound = checkTask.getFormula().getUpperBound().evaluateAsInt(); + STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive."); + if (checkTask.getFormula().isUpperBoundStrict()) { + STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Expected a strict upper step bound that is greater than zero."); + --(*stepBound); + } + STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive."); + + // get the results for the subformulas + storm::modelchecker::SparsePropositionalModelChecker propositionalChecker(this->parametricModel); + STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getLeftSubformula()) && propositionalChecker.canHandle(checkTask.getFormula().getRightSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported"); + storm::storage::BitVector phiStates = std::move(propositionalChecker.check(checkTask.getFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector()); + storm::storage::BitVector psiStates = std::move(propositionalChecker.check(checkTask.getFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector()); + + // get the maybeStates + maybeStates = storm::utility::graph::performProbGreater0(this->parametricModel.getBackwardTransitions(), phiStates, psiStates, true, *stepBound); + maybeStates &= ~psiStates; + + // set the result for all non-maybe states + resultsForNonMaybeStates = std::vector(this->parametricModel.getNumberOfStates(), storm::utility::zero()); + storm::utility::vector::setVectorValues(resultsForNonMaybeStates, psiStates, storm::utility::one()); + + // if there are maybestates, create the parameterLifter + if (!maybeStates.empty()) { + // Create the vector of one-step probabilities to go to target states. + std::vector b = this->parametricModel.getTransitionMatrix().getConstrainedRowSumVector(storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), psiStates); + + parameterLifter = std::make_unique>(this->parametricModel.getTransitionMatrix(), b, maybeStates, maybeStates); + } + + // We know some bounds for the results so set them + lowerResultBound = storm::utility::zero(); + upperResultBound = storm::utility::one(); + } + + template + void SparseDtmcParameterLiftingModelChecker::specifyUntilFormula(CheckTask const& checkTask) { + + // get the results for the subformulas + storm::modelchecker::SparsePropositionalModelChecker propositionalChecker(this->parametricModel); + STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getLeftSubformula()) && propositionalChecker.canHandle(checkTask.getFormula().getRightSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported"); + storm::storage::BitVector phiStates = std::move(propositionalChecker.check(checkTask.getFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector()); + storm::storage::BitVector psiStates = std::move(propositionalChecker.check(checkTask.getFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector()); + + // get the maybeStates + std::pair statesWithProbability01 = storm::utility::graph::performProb01(this->parametricModel.getBackwardTransitions(), phiStates, psiStates); + maybeStates = ~(statesWithProbability01.first | statesWithProbability01.second); + + // set the result for all non-maybe states + resultsForNonMaybeStates = std::vector(this->parametricModel.getNumberOfStates(), storm::utility::zero()); + storm::utility::vector::setVectorValues(resultsForNonMaybeStates, statesWithProbability01.second, storm::utility::one()); + + // if there are maybestates, create the parameterLifter + if (!maybeStates.empty()) { + // Create the vector of one-step probabilities to go to target states. + std::vector b = this->parametricModel.getTransitionMatrix().getConstrainedRowSumVector(storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), psiStates); + + parameterLifter = std::make_unique>(this->parametricModel.getTransitionMatrix(), b, maybeStates, maybeStates); + } + + // We know some bounds for the results so set them + lowerResultBound = storm::utility::zero(); + upperResultBound = storm::utility::one(); + } + + template + void SparseDtmcParameterLiftingModelChecker::specifyReachabilityRewardFormula(CheckTask const& checkTask) { + + // get the results for the subformula + storm::modelchecker::SparsePropositionalModelChecker propositionalChecker(this->parametricModel); + STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported"); + storm::storage::BitVector targetStates = std::move(propositionalChecker.check(checkTask.getFormula().getSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector()); + + // get the maybeStates + storm::storage::BitVector infinityStates = storm::utility::graph::performProb1(this->parametricModel.getBackwardTransitions(), storm::storage::BitVector(this->parametricModel.getNumberOfStates(), true), targetStates); + infinityStates.complement(); + maybeStates = ~(targetStates | infinityStates); + + // set the result for all the non-maybe states + resultsForNonMaybeStates = std::vector(this->parametricModel.getNumberOfStates(), storm::utility::zero()); + storm::utility::vector::setVectorValues(resultsForNonMaybeStates, infinityStates, storm::utility::infinity()); + + // if there are maybestates, create the parameterLifter + if (!maybeStates.empty()) { + // Create the reward vector + STORM_LOG_THROW((checkTask.isRewardModelSet() && this->parametricModel.hasRewardModel(checkTask.getRewardModel())) || (!checkTask.isRewardModelSet() && this->parametricModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model specified by the CheckTask is not available in the given model."); + + typename SparseModelType::RewardModelType const& rewardModel = checkTask.isRewardModelSet() ? this->parametricModel.getRewardModel(checkTask.getRewardModel()) : this->parametricModel.getUniqueRewardModel(); + + std::vector b = rewardModel.getTotalRewardVector(this->parametricModel.getTransitionMatrix()); + + parameterLifter = std::make_unique>(this->parametricModel.getTransitionMatrix(), b, maybeStates, maybeStates); + } + + // We only know a lower bound for the result + lowerResultBound = storm::utility::zero(); + + } + + template + void SparseDtmcParameterLiftingModelChecker::specifyCumulativeRewardFormula(CheckTask const& checkTask) { + + // Obtain the stepBound + stepBound = checkTask.getFormula().getBound().evaluateAsInt(); + if (checkTask.getFormula().isBoundStrict()) { + STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Expected a strict upper step bound that is greater than zero."); + --(*stepBound); + } + STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive."); + + // Create the reward vector + STORM_LOG_THROW((checkTask.isRewardModelSet() && this->parametricModel.hasRewardModel(checkTask.getRewardModel())) || (!checkTask.isRewardModelSet() && this->parametricModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model specified by the CheckTask is not available in the given model."); + typename SparseModelType::RewardModelType const& rewardModel = checkTask.isRewardModelSet() ? this->parametricModel.getRewardModel(checkTask.getRewardModel()) : this->parametricModel.getUniqueRewardModel(); + std::vector b = rewardModel.getTotalRewardVector(this->parametricModel.getTransitionMatrix()); + + parameterLifter = std::make_unique>(this->parametricModel.getTransitionMatrix(), b, storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getColumnCount(), true)); + + + // We only know a lower bound for the result + lowerResultBound = storm::utility::zero(); + } + + template + std::unique_ptr SparseDtmcParameterLiftingModelChecker::computeQuantitativeValues(ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters) { + + parameterLifter->specifyRegion(region, dirForParameters); + + // Set up the solver + auto solver = solverFactory->create(parameterLifter->getMatrix()); + if(lowerResultBound) solver->setLowerBound(lowerResultBound.get()); + if(upperResultBound) solver->setUpperBound(upperResultBound.get()); + if(storm::solver::minimize(dirForParameters) && minSched && !stepBound) solver->setSchedulerHint(std::move(minSched.get())); + if(storm::solver::maximize(dirForParameters) && maxSched && !stepBound) solver->setSchedulerHint(std::move(maxSched.get())); + + // Invoke the solver + if(stepBound) { + assert(*stepBound > 0); + x = std::vector(maybeStates.getNumberOfSetBits(), storm::utility::zero()); + solver->repeatedMultiply(dirForParameters, x, ¶meterLifter->getVector(), *stepBound); + } else { + x.resize(maybeStates.getNumberOfSetBits(), storm::utility::zero()); + solver->solveEquations(dirForParameters, x, parameterLifter->getVector()); + if(storm::solver::minimize(dirForParameters)) { + minSched = std::move(*solver->getScheduler()); + } else { + maxSched = std::move(*solver->getScheduler()); + } + } + + // Get the result for the complete model (including maybestates) + std::vector result = resultsForNonMaybeStates; + auto maybeStateResIt = x.begin(); + for(auto const& maybeState : maybeStates) { + result[maybeState] = *maybeStateResIt; + ++maybeStateResIt; + } + return std::make_unique>(std::move(result)); + } + + + template class SparseDtmcParameterLiftingModelChecker, double>; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.h b/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.h new file mode 100644 index 000000000..d057e570e --- /dev/null +++ b/src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.h @@ -0,0 +1,47 @@ +#pragma once + +#include +#include +#include + +#include "storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h" +#include "storm/storage/BitVector.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" +#include "storm/transformer/ParameterLifter.h" + +namespace storm { + namespace modelchecker { + namespace parametric { + + template + class SparseDtmcParameterLiftingModelChecker : public SparseParameterLiftingModelChecker { + public: + SparseDtmcParameterLiftingModelChecker(SparseModelType const& parametricModel); + SparseDtmcParameterLiftingModelChecker(SparseModelType const& parametricModel, std::unique_ptr>&& solverFactory); + + protected: + + virtual void specifyBoundedUntilFormula(CheckTask const& checkTask) override; + virtual void specifyUntilFormula(CheckTask const& checkTask) override; + virtual void specifyReachabilityRewardFormula(CheckTask const& checkTask) override; + virtual void specifyCumulativeRewardFormula(CheckTask const& checkTask) override; + + virtual std::unique_ptr computeQuantitativeValues(ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters) override; + + + private: + storm::storage::BitVector maybeStates; + std::vector resultsForNonMaybeStates; + boost::optional stepBound; + + std::unique_ptr> parameterLifter; + std::unique_ptr> solverFactory; + + // Results from the most recent solver call. + boost::optional minSched, maxSched; + std::vector x; + boost::optional lowerResultBound, upperResultBound; + }; + } + } +} diff --git a/src/storm/modelchecker/parametric/SparseInstantiationModelChecker.cpp b/src/storm/modelchecker/parametric/SparseInstantiationModelChecker.cpp index 205fea5fe..5fbbf9b6b 100644 --- a/src/storm/modelchecker/parametric/SparseInstantiationModelChecker.cpp +++ b/src/storm/modelchecker/parametric/SparseInstantiationModelChecker.cpp @@ -20,8 +20,7 @@ namespace storm { template void SparseInstantiationModelChecker::specifyFormula(storm::modelchecker::CheckTask const& checkTask) { currentFormula = checkTask.getFormula().asSharedPointer(); - currentCheckTask = std::make_unique>(*currentFormula, checkTask.isOnlyInitialStatesRelevantSet()); - currentCheckTask->setProduceSchedulers(checkTask.isProduceSchedulersSet()); + currentCheckTask = std::make_unique>(checkTask.substituteFormula(*currentFormula).template convertValueType()); } template diff --git a/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.cpp b/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.cpp new file mode 100644 index 000000000..3041a642f --- /dev/null +++ b/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.cpp @@ -0,0 +1,92 @@ +#include "SparseParameterLiftingModelChecker.h" + +#include "storm/adapters/CarlAdapter.h" +#include "storm/logic/FragmentSpecification.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/models/sparse/Dtmc.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" + +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/exceptions/NotSupportedException.h" + +namespace storm { + namespace modelchecker { + namespace parametric { + + template + SparseParameterLiftingModelChecker::SparseParameterLiftingModelChecker(SparseModelType const& parametricModel) : parametricModel(parametricModel) { + //Intentionally left empty + } + + template + void SparseParameterLiftingModelChecker::specifyFormula(storm::modelchecker::CheckTask const& checkTask) { + + currentFormula = checkTask.getFormula().asSharedPointer(); + currentCheckTask = std::make_unique>(checkTask.substituteFormula(*currentFormula).template convertValueType()); + + if(currentCheckTask->getFormula().isProbabilityOperatorFormula()) { + auto const& probOpFormula = currentCheckTask->getFormula().asProbabilityOperatorFormula(); + if(probOpFormula.getSubformula().isBoundedUntilFormula()) { + specifyBoundedUntilFormula(currentCheckTask->substituteFormula(probOpFormula.getSubformula().asBoundedUntilFormula())); + } else if(probOpFormula.getSubformula().isUntilFormula()) { + specifyUntilFormula(currentCheckTask->substituteFormula(probOpFormula.getSubformula().asUntilFormula())); + } else if (probOpFormula.getSubformula().isEventuallyFormula()) { + specifyReachabilityProbabilityFormula(currentCheckTask->substituteFormula(probOpFormula.getSubformula().asEventuallyFormula())); + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property."); + } + } else if (currentCheckTask->getFormula().isRewardOperatorFormula()) { + auto const& rewOpFormula = currentCheckTask->getFormula().asRewardOperatorFormula(); + if(rewOpFormula.getSubformula().isEventuallyFormula()) { + specifyReachabilityRewardFormula(currentCheckTask->substituteFormula(rewOpFormula.getSubformula().asEventuallyFormula())); + } else if (rewOpFormula.getSubformula().isCumulativeRewardFormula()) { + specifyCumulativeRewardFormula(currentCheckTask->substituteFormula(rewOpFormula.getSubformula().asCumulativeRewardFormula())); + } + } + } + + template + std::unique_ptr SparseParameterLiftingModelChecker::check(ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters) { + auto quantitativeResult = computeQuantitativeValues(region, dirForParameters); + if(currentCheckTask->getFormula().hasQuantitativeResult()) { + return quantitativeResult; + } else { + return quantitativeResult->template asExplicitQuantitativeCheckResult().compareAgainstBound(this->currentCheckTask->getFormula().asOperatorFormula().getComparisonType(), this->currentCheckTask->getFormula().asOperatorFormula().template getThresholdAs()); + } + } + + template + void SparseParameterLiftingModelChecker::specifyBoundedUntilFormula(CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property."); + } + + template + void SparseParameterLiftingModelChecker::specifyUntilFormula(CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property."); + } + + template + void SparseParameterLiftingModelChecker::specifyReachabilityProbabilityFormula(CheckTask const& checkTask) { + // transform to until formula + auto untilFormula = std::make_shared(storm::logic::Formula::getTrueFormula(), checkTask.getFormula().getSubformula().asSharedPointer()); + specifyUntilFormula(currentCheckTask->substituteFormula(*untilFormula)); + } + + template + void SparseParameterLiftingModelChecker::specifyReachabilityRewardFormula(CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property."); + } + + template + void SparseParameterLiftingModelChecker::specifyCumulativeRewardFormula(CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property."); + } + + template class SparseParameterLiftingModelChecker, double>; + template class SparseParameterLiftingModelChecker, double>; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h b/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h new file mode 100644 index 000000000..0264a6ccf --- /dev/null +++ b/src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h @@ -0,0 +1,55 @@ +#pragma once + +#include "storm/logic/Formulas.h" +#include "storm/modelchecker/CheckTask.h" +#include "storm/modelchecker/results/CheckResult.h" +#include "storm/modelchecker/parametric/ParameterRegion.h" +#include "storm/solver/OptimizationDirection.h" +#include "storm/utility/parametric.h" + +namespace storm { + namespace modelchecker { + namespace parametric { + + /*! + * Class to approximatively check a formula on a parametric model for all parameter valuations within a region + * It is assumed that all considered valuations are graph-preserving and well defined, i.e., + * * all non-const transition probabilities evaluate to some non-zero value + * * the sum of all outgoing transitions is one + */ + template + class SparseParameterLiftingModelChecker { + public: + SparseParameterLiftingModelChecker(SparseModelType const& parametricModel); + + void specifyFormula(CheckTask const& checkTask); + + /*! + * Checks the specified formula on the given region by applying parameter lifting (Parameter choices are lifted to nondeterministic choices) + * + * @param region the region on which parameter lifting is applied + * @param dirForParameters The optimization direction for the parameter choices. If this is, e.g., minimize, then the returned result will be a lower bound for all results induced by the parameter evaluations inside the region. + */ + std::unique_ptr check(ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters); + + protected: + + virtual void specifyBoundedUntilFormula(CheckTask const& checkTask); + virtual void specifyUntilFormula(CheckTask const& checkTask); + virtual void specifyReachabilityProbabilityFormula(CheckTask const& checkTask); + virtual void specifyReachabilityRewardFormula(CheckTask const& checkTask); + virtual void specifyCumulativeRewardFormula(CheckTask const& checkTask); + + virtual std::unique_ptr computeQuantitativeValues(ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters) = 0; + + + SparseModelType const& parametricModel; + std::unique_ptr> currentCheckTask; + + private: + // store the current formula. Note that currentCheckTask only stores a reference to the formula. + std::shared_ptr currentFormula; + }; + } + } +} diff --git a/src/storm/modelchecker/region/SparseMdpRegionModelChecker.cpp b/src/storm/modelchecker/region/SparseMdpRegionModelChecker.cpp index f5b6d97a0..2a80133ef 100644 --- a/src/storm/modelchecker/region/SparseMdpRegionModelChecker.cpp +++ b/src/storm/modelchecker/region/SparseMdpRegionModelChecker.cpp @@ -319,7 +319,7 @@ namespace storm { template - void SparseRegionModelChecker::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr formula) { + void SparseMdpRegionModelChecker::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr formula) { std::cout << "approx model for mdps not implemented" << std::endl; } diff --git a/src/storm/transformer/ParameterLifter.cpp b/src/storm/transformer/ParameterLifter.cpp new file mode 100644 index 000000000..6ace8a428 --- /dev/null +++ b/src/storm/transformer/ParameterLifter.cpp @@ -0,0 +1,248 @@ +#include "storm/transformer/ParameterLifter.h" + + +#include "storm/adapters/CarlAdapter.h" + +#include "storm/utility/vector.h" + + +namespace storm { + namespace transformer { + + template + ParameterLifter::ParameterLifter(storm::storage::SparseMatrix const& pMatrix, std::vector const& pVector, storm::storage::BitVector const& selectedRows, storm::storage::BitVector const& selectedColumns) { + + // get a mapping from old column indices to new once + std::vector oldToNewColumnIndexMapping(selectedColumns.size(), selectedColumns.size()); + uint_fast64_t newIndex = 0; + for (auto const& oldColumn : selectedColumns) { + oldToNewColumnIndexMapping[oldColumn] = newIndex++; + } + + // Stores which entries of the original matrix/vector are non-constant + storm::storage::BitVector nonConstMatrixEntries(pMatrix.getEntryCount(), false); + storm::storage::BitVector nonConstVectorEntries(pVector.size(), false); + uint_fast64_t pMatrixEntryIndex = 0; + uint_fast64_t pVectorEntryIndex = 0; + + // The matrix builder for the new matrix. The correct number of rows and entries is not known yet. + storm::storage::SparseMatrixBuilder builder(0, pMatrix.getColumnCount(), 0, true, true, pMatrix.getRowCount()); + uint_fast64_t newRowIndex = 0; + for (uint_fast64_t rowIndex = 0; rowIndex < pMatrix.getRowCount(); ++rowIndex) { + builder.newRowGroup(newRowIndex); + + // Gather the occurring variables within this row and set which entries are non-constant + std::set occurringVariables; + for (auto const& entry : pMatrix.getRow(rowIndex)) { + if (!storm::utility::isConstant(entry.getValue())) { + storm::utility::parametric::gatherOccurringVariables(entry.getValue(), occurringVariables); + nonConstMatrixEntries.set(pMatrixEntryIndex++, true); + } + } + + ParametricType const& pVectorEntry = pVector[rowIndex]; + std::set vectorEntryVariables; + if (!storm::utility::isConstant(pVectorEntry)) { + storm::utility::parametric::gatherOccurringVariables(pVectorEntry, vectorEntryVariables); + nonConstVectorEntries.set(pVectorEntryIndex++, true); + } + + // Compute the (abstract) valuation for each row + auto rowValuations = getVerticesOfAbstractRegion(occurringVariables); + + for (auto const& val : rowValuations) { + // Insert matrix entries for each valuation. For non-constant entries, a dummy value is inserted and the function and the valuation are collected. + // The placeholder for the collected function/valuation are stored in the matrixAssignment. The matrixAssignment is completed after the matrix is finished + for (auto const& entry: pMatrix.getRow(rowIndex)) { + if(storm::utility::isConstant(entry.getValue())) { + builder.addNextValue(newRowIndex, oldToNewColumnIndexMapping[entry.getColumn()], storm::utility::convertNumber(entry.getValue())); + } else { + builder.addNextValue(newRowIndex, oldToNewColumnIndexMapping[entry.getColumn()], storm::utility::one()); + ConstantType& placeholder = functionValuationCollector.add(entry.getValue(), val); + matrixAssignment.push_back(std::pair::iterator, ConstantType&>(matrix.begin(), placeholder)); + } + } + + //Insert the vector entry for this row + if(storm::utility::isConstant(pVectorEntry)) { + vector.push_back(storm::utility::convertNumber(pVectorEntry)); + } else { + vector.push_back(storm::utility::one()); + AbstractValuation vectorVal(val); + for(auto const& vectorVar : vectorEntryVariables) { + if(occurringVariables.find(vectorVar) == occurringVariables.end()) { + vectorVal.addParameterUnspecified(vectorVar); + } + } + ConstantType& placeholder = functionValuationCollector.add(pVectorEntry, vectorVal); + vectorAssignment.push_back(std::pair::iterator, ConstantType&>(vector.begin(), placeholder)); + } + ++newRowIndex; + } + } + + // Matrix and vector are now filled with constant results from constant functions and place holders for non-constant functions. + matrix = builder.build(); + vector.shrink_to_fit(); + matrixAssignment.shrink_to_fit(); + vectorAssignment.shrink_to_fit(); + + // Now insert the correct iterators for the matrix and vector assignment + auto matrixAssignmentIt = matrixAssignment.begin(); + uint_fast64_t startEntryOfRow = 0; + for (uint_fast64_t pMatrixRow = 0; pMatrixRow < pMatrix.getRowCount(); ++pMatrixRow) { + uint_fast64_t startEntryOfNextRow = startEntryOfRow + pMatrix.getRow(pMatrixRow).getNumberOfEntries(); + for (uint_fast64_t matrixRow = matrix.getRowGroupIndices()[pMatrixRow]; matrixRow < matrix.getRowGroupIndices()[pMatrixRow + 1]; ++matrixRow) { + auto matrixEntryIt = matrix.getRow(matrixRow).begin(); + for(uint_fast64_t nonConstEntryIndex = nonConstMatrixEntries.getNextSetIndex(startEntryOfRow); nonConstEntryIndex < startEntryOfNextRow; nonConstEntryIndex = nonConstMatrixEntries.getNextSetIndex(nonConstEntryIndex + 1)) { + STORM_LOG_ASSERT(!storm::utility::isConstant((pMatrix.begin() + nonConstEntryIndex)->getValue()), "Expected a non-constant matrix entry."); + matrixAssignmentIt->first = matrixEntryIt + (nonConstEntryIndex - startEntryOfRow); + ++matrixAssignmentIt; + } + } + startEntryOfRow = startEntryOfNextRow; + } + STORM_LOG_ASSERT(matrixAssignmentIt == matrixAssignment.end(), "Unexpected number of entries in the matrix assignment."); + + auto vectorAssignmentIt = vectorAssignment.begin(); + for(auto const& nonConstVectorEntry : nonConstVectorEntries) { + for (uint_fast64_t vectorIndex = matrix.getRowGroupIndices()[nonConstVectorEntry]; vectorIndex != matrix.getRowGroupIndices()[nonConstVectorEntry + 1]; ++vectorIndex) { + vectorAssignmentIt->first = vector.begin() + vectorIndex; + ++vectorAssignmentIt; + } + STORM_LOG_ASSERT(!storm::utility::isConstant(pVector[nonConstVectorEntry]), "Expected a non-constant vector entry."); + } + STORM_LOG_ASSERT(vectorAssignmentIt == vectorAssignment.end(), "Unexpected number of entries in the vector assignment."); + + } + + template + void ParameterLifter::specifyRegion(storm::modelchecker::parametric::ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters) { + // write the evaluation result of each function,evaluation pair into the placeholders + functionValuationCollector.evaluateCollectedFunctions(region, dirForParameters); + + //apply the matrix and vector assignments to write the contents of the placeholder into the matrix/vector + + for(auto& assignment : matrixAssignment) { + assignment.first->setValue(assignment.second); + } + for(auto& assignment : vectorAssignment) { + *assignment.first = assignment.second; + } + } + + template + storm::storage::SparseMatrix const& ParameterLifter::getMatrix() const { + return matrix; + } + + template + std::vector const& ParameterLifter::getVector() const { + return vector; + } + + template + std::vector::AbstractValuation> ParameterLifter::getVerticesOfAbstractRegion(std::set const& variables) const { + std::size_t const numOfVertices = std::pow(2, variables.size()); + std::vector result(numOfVertices); + + for (uint_fast64_t vertexId = 0; vertexId < numOfVertices; ++vertexId) { + //interprete vertexId as a bit sequence + //the consideredVariables.size() least significant bits of vertex will always represent the next vertex + //(00...0 = lower boundaries for all variables, 11...1 = upper boundaries for all variables) + uint_fast64_t variableIndex = 0; + for (auto const& variable : variables) { + if ((vertexId >> variableIndex) % 2 == 0) { + result[vertexId].addParameterLower(variable); + } else { + result[vertexId].addParameterUpper(variable); + } + ++variableIndex; + } + } + return result; + } + + template + bool ParameterLifter::AbstractValuation::operator==(AbstractValuation const& other) const { + return this->lowerPars == other.lowerPars && this->upperPars == other.upperPars && this->unspecifiedPars == other.unspecifiedPars; + } + + template + void ParameterLifter::AbstractValuation::addParameterLower(VariableType const& var) { + lowerPars.insert(var); + } + + template + void ParameterLifter::AbstractValuation::addParameterUpper(VariableType const& var) { + upperPars.insert(var); + } + + template + void ParameterLifter::AbstractValuation::addParameterUnspecified(VariableType const& var) { + unspecifiedPars.insert(var); + } + + template + std::size_t ParameterLifter::AbstractValuation::getHashValue() const { + std::size_t seed = 0; + for (auto const& p : lowerPars) { + carl::hash_add(seed, p); + } + for (auto const& p : upperPars) { + carl::hash_add(seed, p); + } + for (auto const& p : unspecifiedPars) { + carl::hash_add(seed, p); + } + return seed; + } + + template + std::vector> ParameterLifter::AbstractValuation::getConcreteValuations(storm::modelchecker::parametric::ParameterRegion const& region) const { + auto result = region.getVerticesOfRegion(unspecifiedPars); + for(auto& valuation : result) { + for (auto const& lowerPar : lowerPars) { + valuation.insert(std::pair(lowerPar, region.getLowerBoundary(lowerPar))); + } + for (auto const& upperPar : upperPars) { + valuation.insert(std::pair(upperPar, region.getUpperBoundary(upperPar))); + } + } + return result; + } + + template + ConstantType& ParameterLifter::FunctionValuationCollector::add(ParametricType const& function, AbstractValuation const& valuation) { + + // insert the function and the valuation + auto functionsIt = collectedFunctions.insert(std::pair(FunctionValuation(function, valuation), storm::utility::one())).first; + //Note that references to elements of an unordered map remain valid after calling unordered_map::insert. + return functionsIt->second; + } + + template + void ParameterLifter::FunctionValuationCollector::evaluateCollectedFunctions(storm::modelchecker::parametric::ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForUnspecifiedParameters) { + for(auto& collectedFunctionValuationPlaceholder : collectedFunctions) { + ParametricType const& function = collectedFunctionValuationPlaceholder.first.first; + AbstractValuation const& abstrValuation = collectedFunctionValuationPlaceholder.first.second; + ConstantType& placeholder = collectedFunctionValuationPlaceholder.second; + + auto concreteValuations = abstrValuation.getConcreteValuations(region); + auto concreteValuationIt = concreteValuations.begin(); + placeholder = storm::utility::convertNumber(storm::utility::parametric::evaluate(function, *concreteValuationIt)); + for(++concreteValuationIt; concreteValuationIt != concreteValuations.end(); ++concreteValuationIt) { + ConstantType currentResult = storm::utility::convertNumber(storm::utility::parametric::evaluate(function, *concreteValuationIt)); + if(storm::solver::minimize(dirForUnspecifiedParameters)) { + placeholder = std::min(placeholder, currentResult); + } else { + placeholder = std::max(placeholder, currentResult); + } + } + } + } + + + template class ParameterLifter; + } +} diff --git a/src/storm/transformer/ParameterLifter.h b/src/storm/transformer/ParameterLifter.h new file mode 100644 index 000000000..00e15fd05 --- /dev/null +++ b/src/storm/transformer/ParameterLifter.h @@ -0,0 +1,133 @@ +#pragma once + +#include +#include +#include +#include + + +#include "storm/storage/BitVector.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/utility/parametric.h" +#include "storm/modelchecker/parametric/ParameterRegion.h" +#include "storm/solver/OptimizationDirection.h" + +namespace storm { + namespace transformer { + + /*! + * This class lifts parameter choices to nondeterminism: + * For each row in the given matrix that considerd #par parameters, the resulting matrix will have one row group consisting of 2^#par rows. + * When specifying a region, each row within the row group is evaluated w.r.t. one vertex of the region. + * The given vector is handled similarly. + * However, if a vector entry considers a parameter that does not occur in the corresponding matrix row, + * the parameter is directly set such that the vector entry is maximized (or minimized, depending on the specified optimization direction). + * + * @note The row grouping of the original matrix is ignored. + */ + template + class ParameterLifter { + public: + + typedef typename storm::utility::parametric::VariableType::type VariableType; + typedef typename storm::utility::parametric::CoefficientType::type CoefficientType; + + /*! + * Lifts the parameter choices to nondeterminisim. The computation is performed on the submatrix specified by the selected rows and columns + * @param pMatrix the parametric matrix + * @param pVector the parametric vector (the vector size should equal the row count of the matrix) + * @param selectedRows a Bitvector that specifies which rows of the matrix and the vector are considered. + * @param selectedColumns a Bitvector that specifies which columns of the matrix are considered. + */ + ParameterLifter(storm::storage::SparseMatrix const& pMatrix, std::vector const& pVector, storm::storage::BitVector const& selectedRows, storm::storage::BitVector const& selectedColumns); + + void specifyRegion(storm::modelchecker::parametric::ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForParameters); + + // Returns the resulting matrix. Should only be called AFTER specifying a region + storm::storage::SparseMatrix const& getMatrix() const; + + // Returns the resulting vector. Should only be called AFTER specifying a region + std::vector const& getVector() const; + + + private: + /* + * We minimize the number of function evaluations by only calling evaluate() once for each unique pair of function and valuation. + * The result of each evaluation is then written to all positions in the matrix (and the vector) where the corresponding (function,valuation) occurred. + */ + + /* + * During initialization, the actual regions are not known. Hence, we consider abstract valuations, + * where it is only known whether a parameter will be set to either the lower/upper bound of the region or whether this is unspecified + */ + class AbstractValuation { + public: + AbstractValuation() = default; + AbstractValuation(AbstractValuation const& other) = default; + bool operator==(AbstractValuation const& other) const; + + void addParameterLower(VariableType const& var); + void addParameterUpper(VariableType const& var); + void addParameterUnspecified(VariableType const& var); + + std::size_t getHashValue() const; + + /*! + * Returns the concrete valuation(s) (w.r.t. the provided region) represented by this abstract valuation. + * Note that an abstract valuation represents 2^(#unspecified parameters) many concrete valuations. + */ + std::vector> getConcreteValuations(storm::modelchecker::parametric::ParameterRegion const& region) const; + + private: + std::set lowerPars, upperPars, unspecifiedPars; + }; + + /*! + * Collects all occurring pairs of functions and (abstract) valuations. + * We also store a placeholder for the result of each pair. The result is computed and written into the placeholder whenever a region and optimization direction is specified. + */ + class FunctionValuationCollector { + public: + FunctionValuationCollector() = default; + + /*! + * Adds the provided function and valuation. + * Returns a reference to a placeholder in which the evaluation result will be written upon calling evaluateCollectedFunctions) + */ + ConstantType& add(ParametricType const& function, AbstractValuation const& valuation); + + void evaluateCollectedFunctions(storm::modelchecker::parametric::ParameterRegion const& region, storm::solver::OptimizationDirection const& dirForUnspecifiedParameters); + + private: + // Stores a function and a valuation. The valuation is stored as an index of the collectedValuations-vector. + typedef std::pair FunctionValuation; + class FuncValHash{ + public: + std::size_t operator()(FunctionValuation const& fv) const { + std::size_t seed = 0; + carl::hash_add(seed, fv.first); + carl::hash_add(seed, fv.second.getHashValue()); + return seed; + } + }; + + // Stores the collected functions with the valuations together with a placeholder for the result. + std::unordered_map collectedFunctions; + }; + + FunctionValuationCollector functionValuationCollector; + + // Returns the 2^(variables.size()) vertices of the region + std::vector getVerticesOfAbstractRegion(std::set const& variables) const; + + + storm::storage::SparseMatrix matrix; //The resulting matrix; + std::vector::iterator, ConstantType&>> matrixAssignment; // Connection of matrix entries with placeholders + + std::vector vector; //The resulting vector + std::vector::iterator, ConstantType&>> vectorAssignment; // Connection of vector entries with placeholders + + }; + + } +} diff --git a/src/storm/transformer/SparseParametricModelSimplifier.cpp b/src/storm/transformer/SparseParametricModelSimplifier.cpp index ec35bb45e..1c9fef8d4 100644 --- a/src/storm/transformer/SparseParametricModelSimplifier.cpp +++ b/src/storm/transformer/SparseParametricModelSimplifier.cpp @@ -1,8 +1,5 @@ -#include -#include #include "storm/transformer/SparseParametricModelSimplifier.h" - #include "storm/adapters/CarlAdapter.h" #include "storm/models/sparse/Dtmc.h" diff --git a/src/storm/utility/ModelInstantiator.cpp b/src/storm/utility/ModelInstantiator.cpp index c994b129c..3887bfd33 100644 --- a/src/storm/utility/ModelInstantiator.cpp +++ b/src/storm/utility/ModelInstantiator.cpp @@ -1,10 +1,3 @@ -/* - * File: ModelInstantiator.cpp - * Author: Tim Quatmann - * - * Created on February 23, 2016 - */ - #include "storm/utility/ModelInstantiator.h" #include "storm/models/sparse/StandardRewardModel.h" diff --git a/src/storm/utility/ModelInstantiator.h b/src/storm/utility/ModelInstantiator.h index bfb4f687e..9e14f4620 100644 --- a/src/storm/utility/ModelInstantiator.h +++ b/src/storm/utility/ModelInstantiator.h @@ -1,10 +1,3 @@ -/* - * File: ModelInstantiator.h - * Author: Tim Quatmann - * - * Created on February 23, 2016 - */ - #ifndef STORM_UTILITY_MODELINSTANTIATOR_H #define STORM_UTILITY_MODELINSTANTIATOR_H diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index f2d2b98a8..15b05caad 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -393,6 +393,11 @@ namespace storm { RationalFunction convertNumber(double const& number){ return RationalFunction(carl::rationalize(number)); } + + template<> + double convertNumber(RationalFunction const& number){ + return carl::toDouble(number.constantPart()); + } template<> RationalFunction convertNumber(int_fast64_t const& number){ diff --git a/src/storm/utility/parametric.cpp b/src/storm/utility/parametric.cpp index e1fe06a85..86c69addd 100644 --- a/src/storm/utility/parametric.cpp +++ b/src/storm/utility/parametric.cpp @@ -1,10 +1,3 @@ -/* - * File: parametric.cpp - * Author: Tim Quatmann - * - * Created by Tim Quatmann on 08/03/16. - */ - #include #include "storm/utility/parametric.h" @@ -33,6 +26,11 @@ namespace storm { typename CoefficientType::type getConstantPart(storm::RationalFunction const& function){ return function.constantPart(); } + + template<> + void gatherOccurringVariables(storm::RationalFunction const& function, std::set::type>& variableSet){ + function.gatherVariables(variableSet); + } #endif } } diff --git a/src/storm/utility/parametric.h b/src/storm/utility/parametric.h index d07bb9046..7ffdadceb 100644 --- a/src/storm/utility/parametric.h +++ b/src/storm/utility/parametric.h @@ -41,7 +41,13 @@ namespace storm { */ template typename CoefficientType::type getConstantPart(FunctionType const& function); - + + /*! + * Add all variables that occur in the given function to the the given set + */ + template + void gatherOccurringVariables(FunctionType const& function, std::set::type>& variableSet); + } }