From 046afd3804db9ab8612a44694cb23d47c24018ed Mon Sep 17 00:00:00 2001 From: TimQu Date: Sun, 11 Oct 2015 21:09:42 +0200 Subject: [PATCH] Refactored SamplingModel Former-commit-id: b51ed752b4ce5763d0d2303b0bfd75f73c9a3eab --- .../AbstractSparseRegionModelChecker.cpp | 4 +- src/modelchecker/region/SamplingModel.cpp | 328 ++++++++++-------- src/modelchecker/region/SamplingModel.h | 78 +++-- .../region/SparseMdpRegionModelChecker.cpp | 22 +- .../SparseMdpRegionModelCheckerTest.cpp | 4 +- 5 files changed, 254 insertions(+), 182 deletions(-) diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp index f787207a0..f2afacbde 100644 --- a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp +++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp @@ -141,7 +141,7 @@ namespace storm { initializeSamplingModel(*this->getSimpleModel(), this->getSimpleFormula()); std::map emptySubstitution; this->getSamplingModel()->instantiate(emptySubstitution); - this->constantResult = this->getSamplingModel()->computeValues()->template asExplicitQuantitativeCheckResult().getValueVector()[*this->getSamplingModel()->getModel()->getInitialStates().begin()]; + this->constantResult = this->getSamplingModel()->computeInitialStateValue(); } //some more information for statistics... @@ -394,7 +394,7 @@ namespace storm { return this->constantResult.get(); } this->getSamplingModel()->instantiate(point); - return this->getSamplingModel()->computeValues()->template asExplicitQuantitativeCheckResult().getValueVector()[*this->getSamplingModel()->getModel()->getInitialStates().begin()]; + return this->getSamplingModel()->computeInitialStateValue(); } template diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index 8249e3d62..8e3b03107 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -11,143 +11,186 @@ #include "src/models/sparse/Mdp.h" #include "src/models/ModelType.h" #include "models/sparse/StandardRewardModel.h" -#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" -#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" +#include "src/utility/solver.h" +#include "src/solver/LinearEquationSolver.h" +#include "src/solver/MinMaxLinearEquationSolver.h" #include "src/utility/macros.h" #include "src/utility/region.h" #include "src/utility/vector.h" #include "src/exceptions/UnexpectedException.h" #include "src/exceptions/InvalidArgumentException.h" +#include "storage/dd/CuddBdd.h" namespace storm { namespace modelchecker { namespace region { template - SamplingModel::SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr formula) : formula(formula){ - if(this->formula->isProbabilityOperatorFormula()){ + SamplingModel::SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr formula) : solveGoal(storm::logic::isLowerBound(formula->getComparisonType())){ + //First some simple checks and initializations.. + if(formula->isProbabilityOperatorFormula()){ this->computeRewards=false; - } else if(this->formula->isRewardOperatorFormula()){ + } else if(formula->isRewardOperatorFormula()){ this->computeRewards=true; + STORM_LOG_THROW(parametricModel.getType()==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Sampling with rewards is only implemented for Dtmcs"); STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should be unique"); STORM_LOG_THROW(parametricModel.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should have state rewards only"); } else { - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << this->formula << ". Sampling model only supports eventually or reachability reward formulae."); + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << formula << ". Sampling model only supports eventually or reachability reward formulae."); } - //Start with the probabilities - storm::storage::SparseMatrix probabilityMatrix; - std::vector matrixEntryToEvalTableMapping;// This vector will get an entry for every probability matrix entry. - // For the corresponding matrix entry, it stores the corresponding entry in the probability evaluation table. - TableEntry constantEntry; //this value is stored in the matrixEntrytoEvalTableMapping for every constant matrix entry. (also used for rewards later) - initializeProbabilities(parametricModel, probabilityMatrix, matrixEntryToEvalTableMapping, &constantEntry); - - //Now consider rewards - std::unordered_map> rewardModels; - std::vector rewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping - if(this->computeRewards){ - boost::optional> stateRewards; - initializeRewards(parametricModel, stateRewards, rewardEntryToEvalTableMapping, &constantEntry); - rewardModels.insert(std::pair>("", storm::models::sparse::StandardRewardModel(std::move(stateRewards)))); - } - - //Obtain other model ingredients and the sampling model itself - storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); - boost::optional>> noChoiceLabeling; - switch(parametricModel.getType()){ - case storm::models::ModelType::Dtmc: - this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling)); - break; - case storm::models::ModelType::Mdp: - this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling)); - break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Tried to build a sampling model for an unsupported model type"); - } - - //translate the matrixEntryToEvalTableMapping into the actual probability mapping - auto sampleModelEntry = this->model->getTransitionMatrix().begin(); - auto parModelEntry = parametricModel.getTransitionMatrix().begin(); - for(auto tableEntry : matrixEntryToEvalTableMapping){ - STORM_LOG_THROW(sampleModelEntry->getColumn()==parModelEntry->getColumn(), storm::exceptions::UnexpectedException, "The entries of the given parametric model and the constructed sampling model do not match"); - if(tableEntry == &constantEntry){ - sampleModelEntry->setValue(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(parModelEntry->getValue()))); - } else { - this->probabilityMapping.emplace_back(std::make_pair(&(tableEntry->second), sampleModelEntry)); - } - ++sampleModelEntry; - ++parModelEntry; + STORM_LOG_THROW(parametricModel.hasLabel("target"), storm::exceptions::InvalidArgumentException, "The given Model has no \"target\"-statelabel."); + this->targetStates = parametricModel.getStateLabeling().getStates("target"); + STORM_LOG_THROW(parametricModel.hasLabel("sink"), storm::exceptions::InvalidArgumentException, "The given Model has no \"sink\"-statelabel."); + storm::storage::BitVector sinkStates=parametricModel.getStateLabeling().getStates("sink"); + this->maybeStates = ~(this->targetStates | sinkStates); + STORM_LOG_THROW(parametricModel.getInitialStates().getNumberOfSetBits()==1, storm::exceptions::InvalidArgumentException, "The given model has more or less then one initial state"); + storm::storage::sparse::state_type initialState = *parametricModel.getInitialStates().begin(); + STORM_LOG_THROW(maybeStates.get(initialState), storm::exceptions::InvalidArgumentException, "The value in the initial state of the given model is independent of parameters"); + //The (state-)indices in the equation system will be different from the original ones, as the eq-sys only considers maybestates. + //Hence, we use this vector to translate from old indices to new ones. + std::vector newIndices(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates()); //initialize with some illegal index + std::size_t newIndex=0; + for(auto const& index : maybeStates){ + newIndices[index]=newIndex; + ++newIndex; } - //also do this for the rewards - if(this->computeRewards){ - auto sampleModelStateRewardEntry = this->model->getUniqueRewardModel()->second.getStateRewardVector().begin(); - for(auto tableEntry : rewardEntryToEvalTableMapping){ - if(tableEntry != &constantEntry){ - this->stateRewardMapping.emplace_back(std::make_pair(&(tableEntry->second), sampleModelStateRewardEntry)); - } - ++sampleModelStateRewardEntry; - } + + //Now pre-compute the information for the equation system. + initializeProbabilities(parametricModel, newIndices); + if(this->computeRewards){ + initializeRewards(parametricModel, newIndices); } + this->matrixData.assignment.shrink_to_fit(); + this->vectorData.assignment.shrink_to_fit(); + + this->eqSysResult = std::vector(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one() : ConstantType(0.5)); + this->eqSysInitIndex = newIndices[initialState]; + this->solveGoal = storm::solver::SolveGoal(storm::logic::isLowerBound(formula->getComparisonType())); } - + template - void SamplingModel::initializeProbabilities(ParametricSparseModelType const& parametricModel, - storm::storage::SparseMatrix& probabilityMatrix, - std::vector& matrixEntryToEvalTableMapping, - TableEntry* constantEntry) { - // Run through the rows of the original model and obtain the probability evaluation table, a matrix with dummy entries and the mapping between the two. + void SamplingModel::initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector const& newIndices){ + //First run: get a matrix with dummy entries at the new positions and fill the vectorData + ConstantType dummyValue = storm::utility::one(); bool addRowGroups = parametricModel.getTransitionMatrix().hasNontrivialRowGrouping(); - auto curRowGroup = parametricModel.getTransitionMatrix().getRowGroupIndices().begin(); - storm::storage::SparseMatrixBuilder matrixBuilder(parametricModel.getTransitionMatrix().getRowCount(), - parametricModel.getTransitionMatrix().getColumnCount(), + bool isDtmc = (parametricModel.getType()==storm::models::ModelType::Dtmc); //The equation system for DTMCs need the (I-P)-matrix (i.e., we need diagonal entries) + auto numOfMaybeStates=maybeStates.getNumberOfSetBits(); + storm::storage::SparseMatrixBuilder matrixBuilder(addRowGroups ? parametricModel.getTransitionMatrix().getRowCount() : numOfMaybeStates, + numOfMaybeStates, parametricModel.getTransitionMatrix().getEntryCount(), - true, //force dimensions + false, // no force dimensions addRowGroups, - addRowGroups ? parametricModel.getTransitionMatrix().getRowGroupCount() : 0); - matrixEntryToEvalTableMapping.reserve(parametricModel.getTransitionMatrix().getEntryCount()); - std::size_t numOfNonConstEntries=0; - for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - if(addRowGroups && row==*curRowGroup){ - matrixBuilder.newRowGroup(row); - ++curRowGroup; + addRowGroups ? numOfMaybeStates : 0); + this->vectorData.vector.reserve(numOfMaybeStates); //Important! iterators have to remain valid + this->vectorData.assignment.reserve(numOfMaybeStates); + std::size_t curRow = 0; + for (auto oldRowGroup : this->maybeStates){ + if(addRowGroups){ + matrixBuilder.newRowGroup(curRow); + } + for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){ + bool insertedDiagonalEntry = false; + ParametricType targetProbability = storm::utility::region::getNewFunction(storm::utility::zero()); + for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ + if(this->maybeStates.get(oldEntry.getColumn())){ + //check if we need to insert a diagonal entry + if(isDtmc && !insertedDiagonalEntry && newIndices[oldEntry.getColumn()]>=newIndices[oldRowGroup]){ + if(newIndices[oldEntry.getColumn()]>newIndices[oldRowGroup]){ + // There is no diagonal entry in the original matrix. + // Since we later need the matrix (I-P), we already know that the diagonal entry will be one (=1-0) + matrixBuilder.addNextValue(curRow, newIndices[oldRowGroup], storm::utility::one()); + } + insertedDiagonalEntry=true; + } + //Insert dummy entry + matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyValue); + } + else if(!this->computeRewards && this->targetStates.get(oldEntry.getColumn())){ + targetProbability += oldEntry.getValue(); + } + } + if(isDtmc && !insertedDiagonalEntry){ + // There is no diagonal entry in the original matrix. + // Since we later need the matrix (I-P), we already know that the diagonal entry will be one (=1-0) + matrixBuilder.addNextValue(curRow, newIndices[oldRowGroup], storm::utility::one()); + } + if(!this->computeRewards){ + if(storm::utility::isConstant(targetProbability)){ + this->vectorData.vector.push_back(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability))); + } else { + typename std::unordered_map::iterator functionsIt = this->vectorData.functions.insert(FunctionEntry(targetProbability, dummyValue)).first; + this->vectorData.assignment.emplace_back(std::make_pair(this->vectorData.vector.end(), &(functionsIt->second))); + this->vectorData.vector.push_back(dummyValue); + } + } + ++curRow; } - ConstantType dummyEntry=storm::utility::one(); - for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ - if(storm::utility::isConstant(entry.getValue())){ - matrixEntryToEvalTableMapping.emplace_back(constantEntry); - } else { - ++numOfNonConstEntries; - auto evalTableIt = this->probabilityEvaluationTable.insert(TableEntry(entry.getValue(), storm::utility::zero())).first; - matrixEntryToEvalTableMapping.emplace_back(&(*evalTableIt)); + } + this->matrixData.matrix=matrixBuilder.build(); + + //Now run again through both matrices to get the remaining ingredients of the matrixData. + //Note that we need the matrix (I-P) in case of a dtmc. + this->matrixData.assignment.reserve(this->matrixData.matrix.getEntryCount()); + curRow = 0; + for(auto oldRowGroup : this->maybeStates){ + for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){ + auto eqSysMatrixEntry = this->matrixData.matrix.getRow(curRow).begin(); + for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ + if(this->maybeStates.get(oldEntry.getColumn())){ + if(isDtmc && eqSysMatrixEntry->getColumn()==newIndices[oldRowGroup] && eqSysMatrixEntry->getColumn()!=newIndices[oldEntry.getColumn()]){ + //We are at one of the diagonal entries that have been inserted above and for which there is no entry in the original matrix. + //These have already been set to 1 above, so they do not need to be handled here. + ++eqSysMatrixEntry; + } + STORM_LOG_THROW(eqSysMatrixEntry->getColumn()==newIndices[oldEntry.getColumn()], storm::exceptions::UnexpectedException, "old and new entries do not match"); + if(storm::utility::isConstant(oldEntry.getValue())){ + if(isDtmc){ + if(eqSysMatrixEntry->getColumn()==newIndices[oldRowGroup]){ //Diagonal entries get 1-c + eqSysMatrixEntry->setValue(storm::utility::one() - storm::utility::region::convertNumber(storm::utility::region::getConstantPart(oldEntry.getValue()))); + } else { + eqSysMatrixEntry->setValue(storm::utility::zero() - storm::utility::region::convertNumber(storm::utility::region::getConstantPart(oldEntry.getValue()))); + } + } else { + eqSysMatrixEntry->setValue(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(oldEntry.getValue()))); + } + } else { + typename std::unordered_map::iterator functionsIt; + if(isDtmc){ + if(eqSysMatrixEntry->getColumn()==newIndices[oldRowGroup]){ //Diagonal entries get 1-f(x) + functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::one()-oldEntry.getValue(), dummyValue)).first; + } else { + functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::zero()-oldEntry.getValue(), dummyValue)).first; + } + } else { + functionsIt = this->matrixData.functions.insert(FunctionEntry(oldEntry.getValue(), dummyValue)).first; + } + this->matrixData.assignment.emplace_back(std::make_pair(eqSysMatrixEntry, &(functionsIt->second))); + //Note that references to elements of an unordered map remain valid after calling unordered_map::insert. + } + ++eqSysMatrixEntry; + } } - matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry); - dummyEntry=storm::utility::zero(); + ++curRow; } } - this->probabilityMapping.reserve(numOfNonConstEntries); - probabilityMatrix=matrixBuilder.build(); + this->matrixData.matrix.updateNonzeroEntryCount(); } template - void SamplingModel::initializeRewards(ParametricSparseModelType const& parametricModel, - boost::optional>& stateRewards, - std::vector& rewardEntryToEvalTableMapping, - TableEntry* constantEntry) { - // run through the state reward vector of the parametric model. Constant entries can be set directly. Parametric entries are inserted into the table - std::vector stateRewardsAsVector(parametricModel.getNumberOfStates()); - rewardEntryToEvalTableMapping.reserve(parametricModel.getNumberOfStates()); - std::size_t numOfNonConstEntries=0; - for(std::size_t state=0; state::initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices){ + //Run through the state reward vector... + this->vectorData.vector.reserve(this->maybeStates.getNumberOfSetBits()); //Important! iterators have to remain valid + this->vectorData.assignment.reserve(this->maybeStates.getNumberOfSetBits()); + ConstantType dummyValue = storm::utility::one(); + for(auto state : this->maybeStates){ if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])){ - stateRewardsAsVector[state] = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])); - rewardEntryToEvalTableMapping.emplace_back(constantEntry); + this->vectorData.vector.push_back(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state]))); } else { - ++numOfNonConstEntries; - auto evalTableIt = this->probabilityEvaluationTable.insert(TableEntry(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], storm::utility::zero())).first; - rewardEntryToEvalTableMapping.emplace_back(&(*evalTableIt)); + typename std::unordered_map::iterator functionsIt = this->vectorData.functions.insert(FunctionEntry(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], dummyValue)).first; + this->vectorData.assignment.emplace_back(std::make_pair(this->vectorData.vector.end(), &(functionsIt->second))); + this->vectorData.vector.push_back(dummyValue); } } - this->stateRewardMapping.reserve(numOfNonConstEntries); - stateRewards=std::move(stateRewardsAsVector); } template @@ -155,59 +198,56 @@ namespace storm { //Intentionally left empty } - template - std::shared_ptr> const& SamplingModel::getModel() const { - return this->model; - } - template void SamplingModel::instantiate(std::mapconst& point) { - //write entries into evaluation tables - for(auto& tableEntry : this->probabilityEvaluationTable){ - tableEntry.second=storm::utility::region::convertNumber( - storm::utility::region::evaluateFunction(tableEntry.first, point)); + //write results into the placeholders + for(auto& functionResult : this->matrixData.functions){ + functionResult.second=storm::utility::region::convertNumber( + storm::utility::region::evaluateFunction(functionResult.first, point)); } - for(auto& tableEntry : this->rewardEvaluationTable){ - tableEntry.second=storm::utility::region::convertNumber( - storm::utility::region::evaluateFunction(tableEntry.first, point)); + for(auto& functionResult : this->vectorData.functions){ + functionResult.second=storm::utility::region::convertNumber( + storm::utility::region::evaluateFunction(functionResult.first, point)); } - //write the instantiated values to the matrix according to the mappings - for(auto& mappingPair : this->probabilityMapping){ - mappingPair.second->setValue(*(mappingPair.first)); + //write the instantiated values to the matrix and the vector according to the assignment + for(auto& assignment : this->matrixData.assignment){ + assignment.first->setValue(*(assignment.second)); } - if(this->computeRewards){ - for(auto& mappingPair : this->stateRewardMapping){ - *mappingPair.second=*mappingPair.first; - } + for(auto& assignment : this->vectorData.assignment){ + *assignment.first=*assignment.second; } } template - std::unique_ptr SamplingModel::computeValues() { - std::unique_ptr modelChecker; - switch(this->getModel()->getType()){ - case storm::models::ModelType::Dtmc: - modelChecker = std::make_unique>>(*this->model->template as>()); - break; - case storm::models::ModelType::Mdp: - modelChecker = std::make_unique>>(*this->model->template as>()); - break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Tried to build a sampling model for an unsupported model type"); - } - std::unique_ptr resultPtr; - //perform model checking - boost::optional opDir = storm::logic::isLowerBound(this->formula->getComparisonType()) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize; - if(this->computeRewards){ - resultPtr = modelChecker->computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula()); - } - else { - storm::logic::UntilFormula newFormula(storm::logic::Formula::getTrueFormula(), this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula().getSubformula().asSharedPointer()); - resultPtr = modelChecker->computeUntilProbabilities(newFormula, false, opDir); - } - return resultPtr; + std::vector SamplingModel::computeValues() { + invokeSolver(); + std::vector result(this->maybeStates.size()); + storm::utility::vector::setVectorValues(result, this->maybeStates, this->eqSysResult); + storm::utility::vector::setVectorValues(result, this->targetStates, this->computeRewards ? storm::utility::zero() : storm::utility::one()); + storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity() : storm::utility::zero()); + + return result; + } + + template + ConstantType SamplingModel::computeInitialStateValue() { + invokeSolver(); + return this->eqSysResult[this->eqSysInitIndex]; + } + + template<> + void SamplingModel, double>::invokeSolver(){ + std::unique_ptr> solver = storm::utility::solver::LinearEquationSolverFactory().create(this->matrixData.matrix); + solver->solveEquationSystem(this->eqSysResult, this->vectorData.vector); + } + template<> + void SamplingModel, double>::invokeSolver(){ + std::unique_ptr> solver = storm::solver::configureMinMaxLinearEquationSolver(this->solveGoal, storm::utility::solver::MinMaxLinearEquationSolverFactory(), this->matrixData.matrix); + solver->solveEquationSystem(this->eqSysResult, this->vectorData.vector); } + + #ifdef STORM_HAVE_CARL template class SamplingModel>, double>; diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index 244be97d1..3fec3b086 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -12,11 +12,10 @@ #include #include "src/utility/region.h" - #include "src/logic/Formulas.h" -#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "src/models/sparse/Model.h" #include "src/storage/SparseMatrix.h" +#include "src/solver/SolveGoal.h" namespace storm { namespace modelchecker{ @@ -30,15 +29,14 @@ namespace storm { /*! * Creates a sampling model. + * The given model should have the state-labels + * * "target", labeled on states with reachability probability one (reachability reward zero) + * * "sink", labeled on states from which a target state can not be reached. + * The (single) initial state should be disjoint from these states. (otherwise the result would be independent of the parameters, anyway) */ SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr formula); virtual ~SamplingModel(); - /*! - * returns the underlying model - */ - std::shared_ptr> const& getModel() const; - /*! * Instantiates the underlying model according to the given point */ @@ -48,35 +46,55 @@ namespace storm { * Returns the reachability probabilities (or the expected rewards) for every state according to the current instantiation. * Undefined behavior if model has not been instantiated first! */ - std::unique_ptr computeValues(); - + std::vector computeValues(); + + /*! + * Returns the reachability probability (or the expected rewards) of the initial state. + * Undefined behavior if model has not been instantiated first! + */ + ConstantType computeInitialStateValue(); private: - typedef typename std::unordered_map::value_type TableEntry; - - void initializeProbabilities(ParametricSparseModelType const& parametricModel, storm::storage::SparseMatrix& probabilityMatrix, std::vector& matrixEntryToEvalTableMapping, TableEntry* constantEntry); - void initializeRewards(ParametricSparseModelType const& parametricModel, boost::optional>& stateRewards, std::vector& rewardEntryToEvalTableMapping, TableEntry* constantEntry); - - //The model with which we work - std::shared_ptr> model; - //The formula for which we will compute the values - std::shared_ptr formula; + typedef typename std::unordered_map::value_type FunctionEntry; + void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector const& newIndices); + void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices); + void invokeSolver(); + + //Some designated states in the model + storm::storage::BitVector targetStates, maybeStates; + //The last result of the solving the equation system. Also serves as first guess for the next call. + //Note: eqSysResult.size==maybeStates.numberOfSetBits + std::vector eqSysResult; + //The index which represents the result for the initial state in the eqSysResult vector + std::size_t eqSysInitIndex; //A flag that denotes whether we compute probabilities or rewards bool computeRewards; + //The goal we want to accomplish when solving the eq sys. + storm::solver::SolveGoal solveGoal; + - // We store one (unique) entry for every occurring function. - // Whenever a sampling point is given, we can then evaluate the functions - // and store the result to the target value of this map - std::unordered_map probabilityEvaluationTable; - std::unordered_map rewardEvaluationTable; - - //This Vector connects the probability evaluation table with the probability matrix of the model. - //Vector has one entry for every (non-constant) matrix entry. - //pair.first points to an entry in the evaluation table, - //pair.second is an iterator to the corresponding matrix entry - std::vector::iterator>> probabilityMapping; - std::vector::iterator>> stateRewardMapping; + /* The data required for the equation system, i.e., a matrix and a vector. + * + * We use a map to store one (unique) entry for every occurring function. + * The map points to some ConstantType value which serves as placeholder. + * When instantiating the model, the evaluated result of every function is stored in the corresponding placeholder. + * Finally, there is an assignment that connects every non-constant matrix entry + * with a pointer to the value that, on instantiation, needs to be written in that entry. + * + * This way, it is avoided that the same function is evaluated multiple times. + */ + struct MatrixData { + storm::storage::SparseMatrix matrix; //The matrix itself. + std::unordered_map functions; // the occurring functions together with the corresponding placeholders for the result + std::vector::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders + } matrixData; + struct VectorData { + std::vector vector; //The vector itself. + std::unordered_map functions; // the occurring functions together with the corresponding placeholders for the result + std::vector::iterator, ConstantType*>> assignment; // Connection of vector entries with placeholders + } vectorData; + }; } //namespace region diff --git a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp index 8ae59ce99..6018d22d3 100644 --- a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp +++ b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp @@ -96,12 +96,26 @@ namespace storm { } std::cout << "Found that " << stateCounter << " of " << this->getModel().getNumberOfStates() << " states could be eliminated" << std::endl; - - //TODO: Actually eliminate the states... STORM_LOG_WARN("No simplification of the original model (like elimination of constant transitions) is happening. Will just use a copy of the original model"); - simpleModel = std::make_shared(this->getModel()); //Note: an actual copy is technically not necessary.. but we will do it here.. - simpleFormula = this->getSpecifiedFormula(); + + //Get a new labeling + storm::storage::sparse::state_type numStates = this->getModel().getNumberOfStates(); + storm::storage::BitVector sinkStates = ~(maybeStates | targetStates); + storm::models::sparse::StateLabeling labeling(numStates); + storm::storage::BitVector initLabel(this->getModel().getInitialStates()); + labeling.addLabel("init", std::move(initLabel)); + labeling.addLabel("target", std::move(targetStates)); + labeling.addLabel("sink", std::move(sinkStates)); + //Other ingredients + std::unordered_map noRewardModels; + boost::optional>> noChoiceLabeling; + simpleModel = std::make_shared(this->getModel().getTransitionMatrix(), std::move(labeling), std::move(noRewardModels), std::move(noChoiceLabeling)); + + //Get the simplified formula + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + std::shared_ptr eventuallyFormula(new storm::logic::EventuallyFormula(targetFormulaPtr)); + simpleFormula = std::shared_ptr(new storm::logic::ProbabilityOperatorFormula(this->getSpecifiedFormula()->getComparisonType(), this->getSpecifiedFormulaBound(), eventuallyFormula)); } template diff --git a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp index 0ddd75c9a..225d14289 100644 --- a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp @@ -47,8 +47,8 @@ TEST(SparseMdpRegionModelCheckerTest, coin_Prob) { auto allVioRegion=storm::modelchecker::region::ParameterRegion::parseRegion("0.6<=p<=0.75,0.7<=q<=0.72"); EXPECT_NEAR(0.9512773402, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); - EXPECT_NEAR(0.7455934939, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); - EXPECT_NEAR(0.4187982158, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.7455987332, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); + EXPECT_NEAR(0.41880345311, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.01535089684, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.01711494956, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision()); EXPECT_NEAR(0.004422535374, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision());