diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index 658043c63..2fc8da2ce 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/src/modelchecker/region/ApproximationModel.cpp @@ -16,31 +16,32 @@ namespace storm { template - SparseDtmcRegionModelChecker::ApproximationModel::ApproximationModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : mapping(), evaluationTable(), substitutions(), computeRewards(computeRewards){ + 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 matrix with some dummy entries - storm::storage::SparseMatrixBuilder matrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates()); - typename storm::storage::SparseMatrix::index_type approxModelRow=0; - uint_fast64_t numOfNonConstEntries=0; + // (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; for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - matrixBuilder.newRowGroup(approxModelRow); + probabilityMatrixBuilder.newRowGroup(probMatrixRow); // For (1): Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper} - std::set occurringVariables; + std::set occurringProbVariables; for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ - storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringVariables); + storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables); } - std::size_t numOfSubstitutions=1ull< currSubstitution; std::size_t parameterIndex=0; - for(auto const& parameter : occurringVariables){ + for(auto const& parameter : occurringProbVariables){ if((substitutionId>>parameterIndex)%2==0){ currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER)); } @@ -63,53 +64,169 @@ namespace storm { for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ if(!this->parametricTypeComparator.isConstant(entry.getValue())){ distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex)); - ++numOfNonConstEntries; + ++numOfNonConstProbEntries; } - matrixBuilder.addNextValue(approxModelRow, entry.getColumn(), dummyEntry); + probabilityMatrixBuilder.addNextValue(probMatrixRow, entry.getColumn(), dummyEntry); dummyEntry=storm::utility::zero(); } - ++approxModelRow; + ++probMatrixRow; } } + storm::storage::SparseMatrix probabilityMatrix(probabilityMatrixBuilder.build()); + this->probabilityMapping.reserve(numOfNonConstProbEntries); + + std::cout << "2" << std::endl; + // 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); + } + for( auto const& rewVar : occurringRewVariables){ + if(occurringProbVariables.find(rewVar)!=occurringProbVariables.end()){ + makeStateReward=false; + break; + } + } + if(makeStateReward){ + std::map 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)); + } + 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> noStateRewards; - boost::optional> noTransitionRewards; boost::optional>> noChoiceLabeling; - this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(transitionRewards), std::move(noChoiceLabeling)); - //Get the evaluation table. Note that it remains sorted due to the fact that distinctFuncSubs is sorted + //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 mapping - this->mapping.reserve(numOfNonConstEntries); - approxModelRow=0; + //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 (; approxModelRowmodel->getTransitionMatrix().getRowGroupIndices()[row+1]; ++approxModelRow){ - auto appEntry = this->model->getTransitionMatrix().getRow(approxModelRow).begin(); + 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[approxModelRow], storm::utility::zero()); - auto const tableIt= std::find(evaluationTable.begin(), evaluationTable.end(), searchedTuple); - //auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple); - //auto const& tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple); - if(tableIt==evaluationTable.end()){ - std::cout << "did not found tuple in the table: " << parEntry.getValue() << " substitution " << rowSubstitutions[approxModelRow] << " in parametric model at row " << row << " column " << parEntry.getColumn() << " approximation Row " << approxModelRow << std::endl; - } + } 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"); - mapping.emplace_back(std::make_pair(&(std::get<2>(*tableIt)), appEntry)); + 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)); + } + } + } + } + //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); + } + } + } + } } @@ -151,9 +268,57 @@ namespace storm { ); } //write the instantiated values to the matrix according to the mapping - for(auto& mappingPair : this->mapping){ + for(auto& mappingPair : this->probabilityMapping){ mappingPair.second->setValue(*(mappingPair.first)); } + + if(this->computeRewards){ + //Instantiate the substitutions + std::vector> instantiatedRewardSubs(this->rewardSubstitutions.size()); + for(uint_fast64_t substitutionIndex=0; substitutionIndexrewardSubstitutions.size(); ++substitutionIndex){ + for(std::pair const& sub : this->rewardSubstitutions[substitutionIndex]){ + switch(sub.second){ + case TypeOfBound::LOWER: + instantiatedRewardSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first))); + break; + case TypeOfBound::UPPER: + instantiatedRewardSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getUpperBound(sub.first))); + break; + case TypeOfBound::CHOSEOPTIMAL: + //Insert some dummy value + instantiatedRewardSubs[substitutionIndex].insert(std::make_pair(sub.first, storm::utility::zero())); + break; + default: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound when instantiating a reward substitution. Index: " << substitutionIndex << " TypeOfBound: "<< ((int)sub.second)); + } + } + } + //write entries into evaluation table + for(auto& tableEntry : this->rewardEvaluationTable){ + ConstantType minValue=storm::utility::infinity(); + ConstantType maxValue=storm::utility::zero(); + //Iterate over the different combinations of lower and upper bounds and update the min/max values + auto const& vertices=region.getVerticesOfRegion(this->choseOptimalRewardPars[std::get<1>(tableEntry)]); + for(auto const& vertex : vertices){ + //extend the substitution + for(auto const& sub : vertex){ + instantiatedRewardSubs[std::get<1>(tableEntry)][sub.first]=sub.second; + } + ConstantType currValue = storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction( + std::get<0>(tableEntry), + instantiatedRewardSubs[std::get<1>(tableEntry)] + ) + ); + minValue=std::min(minValue, currValue); + maxValue=std::max(maxValue, currValue); + } + std::get<2>(tableEntry)= minValue; + std::get<3>(tableEntry)= maxValue; + } + //Note: the rewards are written to the model as soon as the optimality type is known (i.e. in computeValues) + } + } template @@ -164,6 +329,27 @@ namespace storm { std::unique_ptr resultPtr; if(this->computeRewards){ storm::logic::ReachabilityRewardFormula reachRewFormula(targetFormulaPtr); + //write the reward values into the model + switch(optimalityType){ + case storm::logic::OptimalityType::Minimize: + for(auto& mappingTriple : this->stateRewardMapping){ + *std::get<2>(mappingTriple)=*std::get<0>(mappingTriple); + } + for(auto& mappingTriple : this->transitionRewardMapping){ + std::get<2>(mappingTriple)->setValue(*(std::get<0>(mappingTriple))); + } + break; + case storm::logic::OptimalityType::Maximize: + for(auto& mappingTriple : this->stateRewardMapping){ + *std::get<2>(mappingTriple)=*std::get<1>(mappingTriple); + } + for(auto& mappingTriple : this->transitionRewardMapping){ + std::get<2>(mappingTriple)->setValue(*(std::get<1>(mappingTriple))); + } + break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given optimality Type is not supported."); + } //perform model checking on the mdp resultPtr = modelChecker.computeReachabilityRewards(reachRewFormula, false, optimalityType); } diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h index d9ceb09f8..540cc75a8 100644 --- a/src/modelchecker/region/ApproximationModel.h +++ b/src/modelchecker/region/ApproximationModel.h @@ -50,26 +50,38 @@ namespace storm { private: enum class TypeOfBound { LOWER, - UPPER - }; - + UPPER, + CHOSEOPTIMAL + }; + //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>> mapping; + 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>> 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 //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> rewardEvaluationTable; //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region) std::vector> substitutions; + //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 + std::vector> rewardSubstitutions; + std::vector> choseOptimalRewardPars; //The Model with which we work std::shared_ptr> model; + //A flag that denotes whether we compute probabilities or rewards bool computeRewards; // comparators that can be used to compare constants. diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index 6958f6ac2..bdacf0693 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -71,7 +71,7 @@ namespace storm { 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"); - probabilityMapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); + this->probabilityMapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); } ++samEntry; } @@ -81,7 +81,7 @@ namespace storm { 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"); - stateRewardMapping.emplace_back(std::make_pair(&(tableIt->second), &(this->model->getStateRewardVector()[state]))); + this->stateRewardMapping.emplace_back(std::make_pair(&(tableIt->second), &(this->model->getStateRewardVector()[state]))); } } } diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index f2ffb0c7d..ca719acb3 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -61,6 +61,7 @@ namespace storm { //The model with which we work std::shared_ptr> model; + //A flag that denotes whether we compute probabilities or rewards bool computeRewards; // comparators that can be used to compare constants.