From af505c7e8927fc759c4ca7ebea7b1653e33f15a0 Mon Sep 17 00:00:00 2001 From: TimQu Date: Thu, 3 Sep 2015 20:10:42 +0200 Subject: [PATCH] Faster and more structured initialization of approx and sampling model Former-commit-id: 34c2253a1bae822fb7f2e4a1b646450f9cc72e1e --- .../region/ApproximationModel.cpp | 362 +++++++++--------- src/modelchecker/region/ApproximationModel.h | 18 +- src/modelchecker/region/SamplingModel.cpp | 157 +++++--- src/modelchecker/region/SamplingModel.h | 8 +- .../region/SparseDtmcRegionModelChecker.cpp | 8 +- src/utility/vector.h | 40 ++ 6 files changed, 341 insertions(+), 252 deletions(-) diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index 5ef19802d..4f1ddeb3c 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/src/modelchecker/region/ApproximationModel.cpp @@ -6,10 +6,11 @@ */ #include "src/modelchecker/region/ApproximationModel.h" -#include "src/modelchecker/region/ParameterRegion.h" -#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h" +#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "exceptions/UnexpectedException.h" +#include "src/utility/vector.h" +#include "src/utility/regions.h" +#include "src/exceptions/UnexpectedException.h" namespace storm { namespace modelchecker { @@ -17,31 +18,98 @@ namespace storm { template SparseDtmcRegionModelChecker::ApproximationModel::ApproximationModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : computeRewards(computeRewards){ - // Run through the rows of the original model and obtain - // (1) the different substitutions (this->substitutions) and the substitution used for every row - std::vector rowSubstitutions; - this->substitutions.emplace_back(std::map()); //we want that the empty substitution is always the first one - // (2) the set of distinct pairs - std::set> distinctFuncSubs; - // (3) the MDP Probability matrix with some dummy entries - storm::storage::SparseMatrixBuilder probabilityMatrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates()); - typename storm::storage::SparseMatrix::index_type probMatrixRow=0; - std::size_t numOfNonConstProbEntries=0; + //Start with the probabilities + storm::storage::SparseMatrix probabilityMatrix; + std::vector rowSubstitutions;// the substitution used in every row (required if rewards are computed) + 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 (more precisely: the position in that table). + //We can later transform this mapping into the desired mapping with iterators + const std::size_t constantEntryIndex=std::numeric_limits::max(); //this value is stored in the matrixEntrytoEvalTableMapping for every constant matrix entry. (also used for rewards later) + initializeProbabilities(parametricModel, probabilityMatrix, rowSubstitutions, matrixEntryToEvalTableMapping, constantEntryIndex); + + //Now consider rewards + boost::optional> stateRewards; + boost::optional> transitionRewards; + std::vector stateRewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping + std::vector transitionRewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping + if(this->computeRewards){ + initializeRewards(parametricModel, probabilityMatrix, rowSubstitutions, stateRewards, transitionRewards, stateRewardEntryToEvalTableMapping, transitionRewardEntryToEvalTableMapping, constantEntryIndex); + } + //Obtain other model ingredients and the approximation model itself + storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); + boost::optional>> noChoiceLabeling; + this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(transitionRewards), std::move(noChoiceLabeling)); + + //translate the matrixEntryToEvalTableMapping into the actual probability mapping + typename storm::storage::SparseMatrix::index_type matrixRow=0; + auto tableIndex=matrixEntryToEvalTableMapping.begin(); for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - probabilityMatrixBuilder.newRowGroup(probMatrixRow); - // For (1): Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper} - std::set occurringProbVariables; + for (; matrixRowmodel->getTransitionMatrix().getRowGroupIndices()[row+1]; ++matrixRow){ + auto approxModelEntry = this->model->getTransitionMatrix().getRow(matrixRow).begin(); + for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){ + if(*tableIndex == constantEntryIndex){ + approxModelEntry->setValue(storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parEntry.getValue()))); + } else { + this->probabilityMapping.emplace_back(std::make_pair(&(std::get<2>(this->probabilityEvaluationTable[*tableIndex])), approxModelEntry)); + } + ++approxModelEntry; + ++tableIndex; + } + } + } + if(this->computeRewards){ + //the same for state and transition rewards + auto approxModelStateRewardEntry = this->model->getStateRewardVector().begin(); + for (std::size_t const& tableIndex : stateRewardEntryToEvalTableMapping){ + if(tableIndex != constantEntryIndex){ + this->stateRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(this->rewardEvaluationTable[tableIndex])), &(std::get<3>(this->rewardEvaluationTable[tableIndex])), approxModelStateRewardEntry)); + } + ++approxModelStateRewardEntry; + } + auto approxModelTransitionRewardEntry = this->model->getTransitionRewardMatrix().begin(); + for (std::size_t const& tableIndex : transitionRewardEntryToEvalTableMapping){ + this->transitionRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(this->rewardEvaluationTable[tableIndex])), &(std::get<3>(this->rewardEvaluationTable[tableIndex])), approxModelTransitionRewardEntry)); + ++approxModelTransitionRewardEntry; + } + //Get the sets of reward parameters that map to "CHOSEOPTIMAL". + this->choseOptimalRewardPars = std::vector>(this->rewardSubstitutions.size()); + for(std::size_t index = 0; indexrewardSubstitutions.size(); ++index){ + for(auto const& sub : this->rewardSubstitutions[index]){ + if(sub.second==TypeOfBound::CHOSEOPTIMAL){ + this->choseOptimalRewardPars[index].insert(sub.first); + } + } + } + } + } + + template + void SparseDtmcRegionModelChecker::ApproximationModel::initializeProbabilities(storm::models::sparse::Dtmcconst& parametricModel, + storm::storage::SparseMatrix& probabilityMatrix, + std::vector& rowSubstitutions, + std::vector& matrixEntryToEvalTableMapping, + std::size_t const& constantEntryIndex) { + // Run through the rows of the original model and obtain the different substitutions, the probability evaluation table, + // an MDP probability matrix with some dummy entries, and the mapping between the two + storm::storage::SparseMatrixBuilder matrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates()); + this->probabilitySubstitutions.emplace_back(std::map()); //we want that the empty substitution is always the first one + std::size_t numOfNonConstEntries=0; + typename storm::storage::SparseMatrix::index_type matrixRow=0; + for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ + matrixBuilder.newRowGroup(matrixRow); + // Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper} + std::set occurringVariables; for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ - storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables); + storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringVariables); } - std::size_t numOfSubstitutions=1ull< currSubstitution; std::size_t parameterIndex=0; - for(auto const& parameter : occurringProbVariables){ + for(auto const& parameter : occurringVariables){ if((substitutionId>>parameterIndex)%2==0){ currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER)); } @@ -50,184 +118,119 @@ namespace storm { } ++parameterIndex; } - //Find the current substitution in this->substitutions (add it if we see this substitution the first time) - std::size_t substitutionIndex = std::find(this->substitutions.begin(),this->substitutions.end(), currSubstitution) - this->substitutions.begin(); - if(substitutionIndex==this->substitutions.size()){ - this->substitutions.emplace_back(std::move(currSubstitution)); - } + std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->probabilitySubstitutions, std::move(currSubstitution)); rowSubstitutions.push_back(substitutionIndex); - //Run again through the row and... - //For (2): add pair for the occuring functions and the current substitution - //For (3): add a row with dummy entries for the current substitution - //Note that this is still executed even if no variables occur. However, we will not put any pair into distninctFuncSubs if substitution is empty. + //For every substitution, run again through the row and add an entry in matrixEntryToEvalTableMapping as well as dummy entries in the matrix + //Note that this is still executed once, even if no parameters occur. ConstantType dummyEntry=storm::utility::one(); for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ - if(!this->parametricTypeComparator.isConstant(entry.getValue())){ - distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex)); - ++numOfNonConstProbEntries; + if(this->parametricTypeComparator.isConstant(entry.getValue())){ + matrixEntryToEvalTableMapping.emplace_back(constantEntryIndex); + } else { + ++numOfNonConstEntries; + std::size_t tableIndex=storm::utility::vector::findOrInsert(this->probabilityEvaluationTable, std::tuple(entry.getValue(), substitutionIndex, storm::utility::zero())); + matrixEntryToEvalTableMapping.emplace_back(tableIndex); } - probabilityMatrixBuilder.addNextValue(probMatrixRow, entry.getColumn(), dummyEntry); + matrixBuilder.addNextValue(matrixRow, entry.getColumn(), dummyEntry); dummyEntry=storm::utility::zero(); } - ++probMatrixRow; + ++matrixRow; } } - storm::storage::SparseMatrix probabilityMatrix(probabilityMatrixBuilder.build()); - this->probabilityMapping.reserve(numOfNonConstProbEntries); + this->probabilityMapping.reserve(numOfNonConstEntries); + probabilityMatrix=matrixBuilder.build(); + } + + template + void SparseDtmcRegionModelChecker::ApproximationModel::initializeRewards(storm::models::sparse::Dtmc const& parametricModel, + storm::storage::SparseMatrix const& probabilityMatrix, + std::vector const& rowSubstitutions, + boost::optional >& stateRewards, + boost::optional >& transitionRewards, + std::vector& stateRewardEntryToEvalTableMapping, + std::vector& transitionRewardEntryToEvalTableMapping, + std::size_t const& constantEntryIndex) { + // run through the state reward vector of the parametric model. + // Constant entries can be set directly. + // For Parametric entries we have two cases: + // (1) make state rewards if the reward is independent of parameters that occur in probability functions. + // (2) make transition rewards otherwise. - // Now obtain transition and staterewards (if required) - boost::optional> stateRewards; - boost::optional> transitionRewards; - std::vector stateRewardSubstituions; - std::vector transitionRewardRowSubstitutions; - std::set> distinctRewardFuncSubs; - if(this->computeRewards){ - this->rewardSubstitutions.emplace_back(std::map()); //we want that the empty substitution is always the first one - //stateRewards - std::vector stateRewardsAsVector(parametricModel.getNumberOfStates()); - stateRewardSubstituions = std::vector(parametricModel.getNumberOfStates(),0); //init with empty substitution (=0) - std::size_t numOfNonConstStateRewEntries=0; - //TransitionRewards - storm::storage::SparseMatrixBuilder rewardMatrixBuilder(probabilityMatrix.getRowCount(), probabilityMatrix.getColumnCount(), 0, true, true, probabilityMatrix.getRowGroupCount()); - transitionRewardRowSubstitutions = std::vector(rowSubstitutions.size(),0); //init with empty substitution (=0) - std::size_t numOfNonConstTransitonRewEntries=0; - for(std::size_t state=0; state occurringRewVariables; - std::set occurringProbVariables; - bool makeStateReward=true; //we make state reward whenever no probability parameter occurs in the state reward function - if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ - //constant reward for this state - stateRewardsAsVector[state]=storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state])); - } else { - //reward depends on parameters. Lets find out if also probability parameters occur here. - //If this is the case, the reward depends on the nondeterministic choices and should be given as transition rewards. - stateRewardsAsVector[state]=storm::utility::zero(); - storm::utility::regions::gatherOccurringVariables(parametricModel.getStateRewardVector()[state], occurringRewVariables); - for(auto const& entry : parametricModel.getTransitionMatrix().getRow(state)){ - storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables); + //stateRewards + std::vector stateRewardsAsVector(parametricModel.getNumberOfStates()); + std::size_t numOfNonConstStateRewEntries=0; + //TransitionRewards + storm::storage::SparseMatrixBuilder matrixBuilder(probabilityMatrix.getRowCount(), probabilityMatrix.getColumnCount(), 0, true, true, probabilityMatrix.getRowGroupCount()); + std::size_t numOfNonConstTransitonRewEntries=0; + this->rewardSubstitutions.emplace_back(std::map()); //we want that the empty substitution is always the first one + + for(std::size_t state=0; state occurringRewVariables; + std::set occurringProbVariables; + bool makeStateReward=true; + if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ + stateRewardsAsVector[state]=storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state])); + stateRewardEntryToEvalTableMapping.emplace_back(constantEntryIndex); + } else { + //reward depends on parameters. Lets find out if probability parameters occur here. + //If this is the case, the reward depends on the nondeterministic choices and should be given as transition rewards. + storm::utility::regions::gatherOccurringVariables(parametricModel.getStateRewardVector()[state], occurringRewVariables); + for(auto const& entry : parametricModel.getTransitionMatrix().getRow(state)){ + storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables); + } + for( auto const& rewVar : occurringRewVariables){ + if(occurringProbVariables.find(rewVar)!=occurringProbVariables.end()){ + makeStateReward=false; + break; } - for( auto const& rewVar : occurringRewVariables){ - if(occurringProbVariables.find(rewVar)!=occurringProbVariables.end()){ - makeStateReward=false; - break; - } + } + if(makeStateReward){ + //Get the corresponding substitution and substitutionIndex + std::map rewardSubstitution; + for(auto const& rewardVar : occurringRewVariables){ + rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL)); } - if(makeStateReward){ + std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->rewardSubstitutions, std::move(rewardSubstitution)); + //insert table entry and dummy data in the stateRewardVector + std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::tuple(parametricModel.getStateRewardVector()[state], substitutionIndex, storm::utility::zero(), storm::utility::zero())); + stateRewardEntryToEvalTableMapping.emplace_back(tableIndex); + stateRewardsAsVector[state]=storm::utility::zero(); + ++numOfNonConstStateRewEntries; + } else { + for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow rewardSubstitution; for(auto const& rewardVar : occurringRewVariables){ - rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL)); - } - //Find the substitution in this->rewardSubstitutions (add it if we see this substitution the first time) - std::size_t substitutionIndex = std::find(this->rewardSubstitutions.begin(),this->rewardSubstitutions.end(), rewardSubstitution) - this->rewardSubstitutions.begin(); - if(substitutionIndex==this->rewardSubstitutions.size()){ - this->rewardSubstitutions.emplace_back(std::move(rewardSubstitution)); - } - distinctRewardFuncSubs.insert(std::make_pair(parametricModel.getStateRewardVector()[state], substitutionIndex)); - ++numOfNonConstStateRewEntries; - stateRewardSubstituions[state]=substitutionIndex; - } else { - for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow()); - } - //Get the corresponding substitution - std::map rewardSubstitution; - for(auto const& rewardVar : occurringRewVariables){ - typename std::map::iterator const substitutionIt=this->substitutions[rowSubstitutions[matrixRow]].find(rewardVar); - if(substitutionIt==this->substitutions[rowSubstitutions[matrixRow]].end()){ - rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL)); - } else { - rewardSubstitution.insert(*substitutionIt); - } - } - //Find the substitution in this->rewardSubstitutions (add it if we see this substitution the first time) - std::size_t substitutionIndex = std::find(this->rewardSubstitutions.begin(),this->rewardSubstitutions.end(), rewardSubstitution) - this->rewardSubstitutions.begin(); - if(substitutionIndex==this->rewardSubstitutions.size()){ - this->rewardSubstitutions.emplace_back(std::move(rewardSubstitution)); + typename std::map::iterator const substitutionIt=this->probabilitySubstitutions[rowSubstitutions[matrixRow]].find(rewardVar); + if(substitutionIt==this->probabilitySubstitutions[rowSubstitutions[matrixRow]].end()){ + rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL)); + } else { + rewardSubstitution.insert(*substitutionIt); } - distinctRewardFuncSubs.insert(std::make_pair(parametricModel.getStateRewardVector()[state], substitutionIndex)); - ++numOfNonConstTransitonRewEntries; - transitionRewardRowSubstitutions[matrixRow]=substitutionIndex; } - } - } - } - stateRewards=std::move(stateRewardsAsVector); - this->stateRewardMapping.reserve(numOfNonConstStateRewEntries); - transitionRewards=rewardMatrixBuilder.build(); - this->transitionRewardMapping.reserve(numOfNonConstTransitonRewEntries); - } - //Obtain other model ingredients and the approximation model itself - storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); - boost::optional>> noChoiceLabeling; - this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(transitionRewards), std::move(noChoiceLabeling)); - - //Get the (probability) evaluation table. Note that it remains sorted due to the fact that distinctFuncSubs is sorted - this->evaluationTable.reserve(distinctFuncSubs.size()); - for(auto const& funcSub : distinctFuncSubs){ - this->evaluationTable.emplace_back(funcSub.first, funcSub.second, storm::utility::zero()); - } - //Fill in the entries for the probability mapping - probMatrixRow=0; - for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - for (; probMatrixRowmodel->getTransitionMatrix().getRowGroupIndices()[row+1]; ++probMatrixRow){ - auto appEntry = this->model->getTransitionMatrix().getRow(probMatrixRow).begin(); - for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){ - if(this->parametricTypeComparator.isConstant(parEntry.getValue())){ - appEntry->setValue(storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parEntry.getValue()))); - } else { - std::tuple searchedTuple(parEntry.getValue(), rowSubstitutions[probMatrixRow], storm::utility::zero()); - auto const tableIt= std::find(this->evaluationTable.begin(), this->evaluationTable.end(), searchedTuple); - STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the evaluationTable. Either the table is missing that tuple or it is not sorted properly"); - this->probabilityMapping.emplace_back(std::make_pair(&(std::get<2>(*tableIt)), appEntry)); - } - ++appEntry; - } - } - } - if(this->computeRewards){ - //Get the (reward) evaluation table. Note that it remains sorted due to the fact that distinctRewFuncSubs is sorted - this->rewardEvaluationTable.reserve(distinctRewardFuncSubs.size()); - for(auto const& funcSub : distinctRewardFuncSubs){ - this->rewardEvaluationTable.emplace_back(funcSub.first, funcSub.second, storm::utility::zero(), storm::utility::zero()); - } - - for (std::size_t state=0; state< parametricModel.getNumberOfStates(); ++state){ - //check if the state reward is not constant - if(stateRewardSubstituions[state]!=0){ - std::tuple searchedTuple(parametricModel.getStateRewardVector()[state], stateRewardSubstituions[state], storm::utility::zero(), storm::utility::zero()); - auto const tableIt= std::find(this->rewardEvaluationTable.begin(), this->rewardEvaluationTable.end(), searchedTuple); - STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the rewardEvaluationTable. Either the table is missing that tuple or it is not sorted properly (stateRewards)"); - this->stateRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(*tableIt)), &(std::get<3>(*tableIt)), &(this->model->getStateRewardVector()[state]))); - } - //check if transitionrewards are not constant - for(auto matrixRow=this->model->getTransitionRewardMatrix().getRowGroupIndices()[state]; matrixRowmodel->getTransitionRewardMatrix().getRowGroupIndices()[state+1]; ++matrixRow){ - if(transitionRewardRowSubstitutions[matrixRow]!=0){ - //Note that all transitions of the same row will have the same reward value - std::tuple searchedTuple(parametricModel.getStateRewardVector()[state], transitionRewardRowSubstitutions[matrixRow], storm::utility::zero(), storm::utility::zero()); - auto const tableIt= std::find(this->rewardEvaluationTable.begin(), this->rewardEvaluationTable.end(), searchedTuple); - STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the rewardEvaluationTable. Either the table is missing that tuple or it is not sorted properly (transitionRewards)"); - for(auto transitionMatrixEntry=this->model->getTransitionRewardMatrix().getRow(matrixRow).begin(); transitionMatrixEntrymodel->getTransitionRewardMatrix().getRow(matrixRow).end(); ++transitionMatrixEntry){ - this->transitionRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(*tableIt)), &(std::get<3>(*tableIt)), transitionMatrixEntry)); + std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->rewardSubstitutions, std::move(rewardSubstitution)); + //insert table entries and dummy data + std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::tuple(parametricModel.getStateRewardVector()[state], substitutionIndex, storm::utility::zero(), storm::utility::zero())); + for(auto const& matrixEntry : probabilityMatrix.getRow(matrixRow)){ + transitionRewardEntryToEvalTableMapping.emplace_back(tableIndex); + matrixBuilder.addNextValue(matrixRow, matrixEntry.getColumn(), storm::utility::zero()); + ++numOfNonConstTransitonRewEntries; } } - } - } - //Get the sets of reward parameters that map to "CHOSEOPTIMAL". - this->choseOptimalRewardPars = std::vector>(this->rewardSubstitutions.size()); - for(std::size_t index = 0; indexrewardSubstitutions.size(); ++index){ - for(auto const& sub : this->rewardSubstitutions[index]){ - if(sub.second==TypeOfBound::CHOSEOPTIMAL){ - this->choseOptimalRewardPars[index].insert(sub.first); - } + stateRewardsAsVector[state]=storm::utility::zero(); + stateRewardEntryToEvalTableMapping.emplace_back(constantEntryIndex); } } } + stateRewards=std::move(stateRewardsAsVector); + this->stateRewardMapping.reserve(numOfNonConstStateRewEntries); + transitionRewards=matrixBuilder.build(); + this->transitionRewardMapping.reserve(numOfNonConstTransitonRewEntries); } + template SparseDtmcRegionModelChecker::ApproximationModel::~ApproximationModel() { //Intentionally left empty @@ -241,9 +244,9 @@ namespace storm { template void SparseDtmcRegionModelChecker::ApproximationModel::instantiate(const ParameterRegion& region) { //Instantiate the substitutions - std::vector> instantiatedSubs(this->substitutions.size()); - for(uint_fast64_t substitutionIndex=0; substitutionIndexsubstitutions.size(); ++substitutionIndex){ - for(std::pair const& sub : this->substitutions[substitutionIndex]){ + std::vector> instantiatedSubs(this->probabilitySubstitutions.size()); + for(uint_fast64_t substitutionIndex=0; substitutionIndexprobabilitySubstitutions.size(); ++substitutionIndex){ + for(std::pair const& sub : this->probabilitySubstitutions[substitutionIndex]){ switch(sub.second){ case TypeOfBound::LOWER: instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first))); @@ -257,7 +260,7 @@ namespace storm { } } //write entries into evaluation table - for(auto& tableEntry : this->evaluationTable){ + for(auto& tableEntry : this->probabilityEvaluationTable){ std::get<2>(tableEntry)=storm::utility::regions::convertNumber( storm::utility::regions::evaluateFunction( std::get<0>(tableEntry), @@ -316,7 +319,6 @@ namespace storm { } //Note: the rewards are written to the model as soon as the optimality type is known (i.e. in computeValues) } - } template diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h index 540cc75a8..e7037ddae 100644 --- a/src/modelchecker/region/ApproximationModel.h +++ b/src/modelchecker/region/ApproximationModel.h @@ -10,6 +10,8 @@ #include "src/modelchecker/region/SparseDtmcRegionModelChecker.h" #include "src/models/sparse/Mdp.h" +#include "src/modelchecker/region/ParameterRegion.h" +#include "src/storage/SparseMatrix.h" namespace storm { namespace modelchecker { @@ -54,24 +56,28 @@ namespace storm { CHOSEOPTIMAL }; + void initializeProbabilities(storm::models::sparse::Dtmc const& parametricModel, storm::storage::SparseMatrix& probabilityMatrix, std::vector& rowSubstitutions, std::vector& matrixEntryToEvalTableMapping, std::size_t const& constantEntryIndex); + void initializeRewards(storm::models::sparse::Dtmc const& parametricModel, storm::storage::SparseMatrix const& probabilityMatrix, std::vector const& rowSubstitutions, boost::optional>& stateRewards, boost::optional>& transitionRewards, std::vector& stateRewardEntryToEvalTableMapping, std::vector& transitionRewardEntryToEvalTableMapping, std::size_t const& constantEntryIndex); + //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; //similar for the rewards. But now the first entry points to a minimal and the second one to a maximal value. //The third entry points to the state reward vector and the transitionRewardMatrix, respectively. - std::vector> stateRewardMapping; + std::vector::iterator>> stateRewardMapping; std::vector::iterator>> transitionRewardMapping; - //Vector has one entry for - //(every distinct, non-constant function that occurs somewhere in the model) x (the required combinations of lower and upper bounds of the region) - //The second entry represents a substitution as an index in the substitutions vector + //Vector has one (unique) entry for every occurring pair of a non-constant function and + // a substitution, i.e., a mapping of variables to a TypeOfBound + //The second entry represents the substitution as an index in the substitutions vector //The third entry should contain the result when evaluating the function in the first entry, regarding the substitution given by the second entry. - std::vector> evaluationTable; + std::vector> probabilityEvaluationTable; + //For rewards, we have the third entry for the minimal value and the fourth entry for the maximal value std::vector> rewardEvaluationTable; //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region) - std::vector> substitutions; + std::vector> probabilitySubstitutions; //Same for the different substitutions for the reward functions. //In addition, we store the parameters for which the correct substitution //depends on the region and whether to minimize/maximize diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index bdacf0693..370070392 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -6,85 +6,113 @@ */ #include "src/modelchecker/region/SamplingModel.h" -#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h" +#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "exceptions/UnexpectedException.h" +#include "src/utility/vector.h" +#include "src/utility/regions.h" +#include "src/exceptions/UnexpectedException.h" namespace storm { namespace modelchecker { template - SparseDtmcRegionModelChecker::SamplingModel::SamplingModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : probabilityMapping(), stateRewardMapping(), evaluationTable(), computeRewards(computeRewards){ - // Run through the rows of the original model and obtain the set of distinct functions as well as a matrix with dummy entries - std::set functionSet; - storm::storage::SparseMatrixBuilder matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount()); - std::size_t numOfNonConstProbEntries=0; - for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - ConstantType dummyEntry=storm::utility::one(); - for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ - if(!this->parametricTypeComparator.isConstant(entry.getValue())){ - functionSet.insert(entry.getValue()); - ++numOfNonConstProbEntries; - } - matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry); - dummyEntry=storm::utility::zero(); - } - } - this->probabilityMapping.reserve(numOfNonConstProbEntries); + SparseDtmcRegionModelChecker::SamplingModel::SamplingModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : probabilityMapping(), stateRewardMapping(), probabilityEvaluationTable(), computeRewards(computeRewards){ + //Start with the probabilities + storm::storage::SparseMatrix probabilityMatrix; + // This vector will get an entry for every probability matrix entry. + // For the corresponding matrix entry, it stores the right entry in the probability evaluation table (more precisely: the position in that table). + std::vector matrixEntryToEvalTableMapping; + const std::size_t constantEntryIndex=std::numeric_limits::max(); //this value is stored in the matrixEntrytoEvalTableMapping for every constant matrix entry. (also used for rewards later) + initializeProbabilities(parametricModel, probabilityMatrix, matrixEntryToEvalTableMapping, constantEntryIndex); - //Obtain other model ingredients and the sampling model itself - storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); + //Now consider rewards boost::optional> stateRewards; + std::vector rewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping if(this->computeRewards){ - std::size_t numOfNonConstRewEntries=0; - std::vector stateRewardsAsVector(parametricModel.getNumberOfStates()); - for(std::size_t state=0; stateparametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ - stateRewardsAsVector[state] = storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state])); - } else { - stateRewardsAsVector[state] = storm::utility::zero(); - functionSet.insert(parametricModel.getStateRewardVector()[state]); - ++numOfNonConstRewEntries; - } - } - stateRewards=std::move(stateRewardsAsVector); - this->stateRewardMapping.reserve(numOfNonConstRewEntries); + initializeRewards(parametricModel, stateRewards, rewardEntryToEvalTableMapping, constantEntryIndex); } + + //Obtain other model ingredients and the sampling model itself + storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); boost::optional> noTransitionRewards; boost::optional>> noChoiceLabeling; - this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); - - //Get the evaluation table. Note that it remains sorted due to the fact that probFunctionSet is sorted - this->evaluationTable.reserve(functionSet.size()); - for(auto const& func : functionSet){ - this->evaluationTable.emplace_back(func, storm::utility::zero()); - } - - //Fill in the entries for the mappings - auto samEntry = this->model->getTransitionMatrix().begin(); - for(auto const& parEntry : parametricModel.getTransitionMatrix()){ - if(this->parametricTypeComparator.isConstant(parEntry.getValue())){ - samEntry->setValue(storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parEntry.getValue()))); + this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + + //translate the matrixEntryToEvalTableMapping into the actual probability mapping + auto sampleModelEntry = this->model->getTransitionMatrix().begin(); + auto parModelEntry = parametricModel.getTransitionMatrix().begin(); + for(std::size_t const& tableIndex : 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(tableIndex == constantEntryIndex){ + sampleModelEntry->setValue(storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parModelEntry->getValue()))); + } else { + this->probabilityMapping.emplace_back(std::make_pair(&(this->probabilityEvaluationTable[tableIndex].second), sampleModelEntry)); } - else { - std::pair searchedPair(parEntry.getValue(), storm::utility::zero()); - auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair); - STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly"); - this->probabilityMapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); + ++sampleModelEntry; + ++parModelEntry; + } + //also do this for the rewards thing + if(this->computeRewards){ + auto sampleModelStateRewardEntry = this->model->getStateRewardVector().begin(); + for(std::size_t const& tableIndex : rewardEntryToEvalTableMapping){ + if(tableIndex != constantEntryIndex){ + this->stateRewardMapping.emplace_back(std::make_pair(&(this->rewardEvaluationTable[tableIndex].second), sampleModelStateRewardEntry)); + } + ++sampleModelStateRewardEntry; } - ++samEntry; } - if(this->computeRewards){ - for(std::size_t state=0; stateparametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ - std::pair searchedPair(parametricModel.getStateRewardVector()[state], storm::utility::zero()); - auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair); - STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly"); - this->stateRewardMapping.emplace_back(std::make_pair(&(tableIt->second), &(this->model->getStateRewardVector()[state]))); + } + + template + void SparseDtmcRegionModelChecker::SamplingModel::initializeProbabilities(storm::models::sparse::Dtmc const& parametricModel, + storm::storage::SparseMatrix& probabilityMatrix, + std::vector& matrixEntryToEvalTableMapping, + std::size_t const& constantEntryIndex) { + // 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. + storm::storage::SparseMatrixBuilder matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount()); + matrixEntryToEvalTableMapping.reserve(parametricModel.getTransitionMatrix().getEntryCount()); + std::size_t numOfNonConstEntries=0; + for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ + ConstantType dummyEntry=storm::utility::one(); + for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ + if(this->parametricTypeComparator.isConstant(entry.getValue())){ + matrixEntryToEvalTableMapping.emplace_back(constantEntryIndex); + } else { + ++numOfNonConstEntries; + std::size_t tableIndex=storm::utility::vector::findOrInsert(this->probabilityEvaluationTable, std::pair(entry.getValue(), storm::utility::zero())); + matrixEntryToEvalTableMapping.emplace_back(tableIndex); } + matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry); + dummyEntry=storm::utility::zero(); + } + } + this->probabilityMapping.reserve(numOfNonConstEntries); + probabilityMatrix=matrixBuilder.build(); + } + + template + void SparseDtmcRegionModelChecker::SamplingModel::initializeRewards(storm::models::sparse::Dtmc const& parametricModel, + boost::optional>& stateRewards, + std::vector& rewardEntryToEvalTableMapping, + std::size_t const& constantEntryIndex) { + // 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; stateparametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ + stateRewardsAsVector[state] = storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state])); + rewardEntryToEvalTableMapping.emplace_back(constantEntryIndex); + } else { + ++numOfNonConstEntries; + stateRewardsAsVector[state] = storm::utility::zero(); + std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::pair(parametricModel.getStateRewardVector()[state], storm::utility::zero())); + rewardEntryToEvalTableMapping.emplace_back(tableIndex); } } + this->stateRewardMapping.reserve(numOfNonConstEntries); + stateRewards=std::move(stateRewardsAsVector); } template @@ -99,11 +127,16 @@ namespace storm { template void SparseDtmcRegionModelChecker::SamplingModel::instantiate(std::mapconst& point) { - //write entries into evaluation table - for(auto& tableEntry : this->evaluationTable){ + //write entries into evaluation tables + for(auto& tableEntry : this->probabilityEvaluationTable){ tableEntry.second=storm::utility::regions::convertNumber( storm::utility::regions::evaluateFunction(tableEntry.first, point)); } + for(auto& tableEntry : this->rewardEvaluationTable){ + tableEntry.second=storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction(tableEntry.first, point)); + } + //write the instantiated values to the matrix according to the mappings for(auto& mappingPair : this->probabilityMapping){ mappingPair.second->setValue(*(mappingPair.first)); diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index ca719acb3..9eb1c2210 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -10,6 +10,7 @@ #include "src/modelchecker/region/SparseDtmcRegionModelChecker.h" #include "src/models/sparse/Dtmc.h" +#include "src/storage/SparseMatrix.h" namespace storm { namespace modelchecker{ @@ -47,16 +48,19 @@ namespace storm { private: + void initializeProbabilities(storm::models::sparse::Dtmc const& parametricModel, storm::storage::SparseMatrix& probabilityMatrix, std::vector& matrixEntryToEvalTableMapping, std::size_t const& constantEntryIndex); + void initializeRewards(storm::models::sparse::Dtmc const& parametricModel, boost::optional>& stateRewards, std::vector& rewardEntryToEvalTableMapping, std::size_t const& constantEntryIndex); //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> stateRewardMapping; + std::vector::iterator>> stateRewardMapping; //Vector has one entry for every distinct, non-constant function that occurs somewhere in the model. //The second entry should contain the result when evaluating the function in the first entry. - std::vector> evaluationTable; + std::vector> probabilityEvaluationTable; + std::vector> rewardEvaluationTable; //The model with which we work std::shared_ptr> model; diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp index d7af4d61e..5bc5ba4a7 100644 --- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp +++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp @@ -91,7 +91,8 @@ namespace storm { initializeApproximationModel(*this->simplifiedModel); } //now create the model used for Sampling - if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE || (storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){ + if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE || + (!storm::settings::regionSettings().doSample() && storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){ initializeSamplingModel(*this->simplifiedModel); } //Check if the reachability function needs to be computed @@ -295,7 +296,7 @@ namespace storm { labeling.addLabel("sink", std::move(sinkLabel)); // other ingredients if(this->computeRewards){ - storm::utility::vector::selectVectorValues(stateRewards.get(), subsystem, stateRewards.get()); + storm::utility::vector::selectVectorValues(stateRewards.get(), subsystem); stateRewards->push_back(storm::utility::zero()); //target state stateRewards->push_back(storm::utility::zero()); //sink state } @@ -379,6 +380,9 @@ namespace storm { // first. storm::utility::vector::selectVectorValues(stateRewards, maybeStates, this->model.getStateRewardVector()); } + for(auto& stateReward: stateRewards){ + storm::utility::simplify(stateReward); + } } diff --git a/src/utility/vector.h b/src/utility/vector.h index 634d2bc5a..d9163d5c8 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -21,6 +21,27 @@ namespace storm { namespace utility { namespace vector { + /*! + * Finds the given element in the given vector. + * If the vector does not contain the element, it is inserted (at the end of vector). + * Either way, the returned value will be the position of the element inside the vector, + * + * @note old indices to other elements remain valid, as the vector will not be sorted. + * + * @param vector The vector in which the element is searched and possibly insert + * @param element The element that will be searched for (or inserted) + * + * @return The position at which the element is located + */ + template + std::size_t findOrInsert(std::vector& vector, T&& element){ + std::size_t position=std::find(vector.begin(), vector.end(), element) - vector.begin(); + if(position==vector.size()){ + vector.emplace_back(std::move(element)); + } + return position; + } + /*! * Sets the provided values at the provided positions in the given vector. * @@ -64,6 +85,25 @@ namespace storm { } } + /*! + * Selects the elements from the given vector at the specified positions and writes them consecutively into the same vector. + * The size of the vector is reduced accordingly. + * @param vector The vector from which to select the elements and into which the selected elements are to be written. + * @param positions The positions at which to select the elements from the values vector. + */ + template + void selectVectorValues(std::vector& vector, storm::storage::BitVector const& positions) { + uint_fast64_t oldPosition = 0; + for (auto position : positions) { + if(oldPosition!=position){ + vector[oldPosition++] = std::move(vector[position]); + } else { + ++oldPosition; + } + } + vector.resize(oldPosition); + } + /*! * Selects groups of elements from a vector at the specified positions and writes them consecutively into another vector. *