Browse Source

refactored approximation model (almost done)

Former-commit-id: c7f285906b
tempestpy_adaptions
TimQu 9 years ago
parent
commit
67ff27954e
  1. 467
      src/modelchecker/region/ApproximationModel.cpp
  2. 106
      src/modelchecker/region/ApproximationModel.h
  3. 106
      src/modelchecker/region/SamplingModel.cpp
  4. 7
      src/modelchecker/region/SamplingModel.h

467
src/modelchecker/region/ApproximationModel.cpp

@ -10,10 +10,11 @@
#include "src/models/sparse/Mdp.h"
#include "src/models/ModelType.h"
#include "src/models/sparse/StandardRewardModel.h"
#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "src/solver/MinMaxLinearEquationSolver.h"
#include "src/solver/GameSolver.h"
#include "src/utility/macros.h"
#include "src/utility/region.h"
#include "src/utility/solver.h"
#include "src/utility/vector.h"
#include "src/exceptions/UnexpectedException.h"
#include "src/exceptions/InvalidArgumentException.h"
@ -24,135 +25,86 @@ namespace storm {
namespace region {
template<typename ParametricSparseModelType, typename ConstantType>
ApproximationModel<ParametricSparseModelType, ConstantType>::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula) : formula(formula){
if(this->formula->isProbabilityOperatorFormula()){
ApproximationModel<ParametricSparseModelType, ConstantType>::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula) : player1Goal(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, "Approximation with rewards is only implemented for Dtmcs");
STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should be unique");
STORM_LOG_THROW(parametricModel.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should have state rewards only");
} else {
STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << this->formula << ". Approximation model only supports eventually or reachability reward formulae.");
}
//Start with the probabilities
storm::storage::SparseMatrix<ConstantType> probabilityMatrix;
std::vector<typename storm::storage::SparseMatrix<ConstantType>::index_type> approxRowGroupIndices;// Indicates which rows of the probabilityMatrix belong to a given row of the parametricModel
std::vector<std::size_t> rowSubstitutions;// the substitution used in every row (required if rewards are computed)
std::vector<ProbTableEntry*> 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.
//We can later transform this mapping into the desired mapping with iterators
ProbTableEntry constantProbEntry; //this value is stored in the matrixEntrytoEvalTableMapping for every constant probability matrix entry.
initializeProbabilities(parametricModel, probabilityMatrix, approxRowGroupIndices, rowSubstitutions, matrixEntryToEvalTableMapping, &constantProbEntry);
//Now consider rewards
std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ConstantType>> rewardModels;
std::vector<RewTableEntry*> rewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping
RewTableEntry constantRewEntry; //this value is stored in the rewardEntryToEvalTableMapping for every constant reward vector entry.
if(this->computeRewards){
std::vector<ConstantType> stateActionRewards;
initializeRewards(parametricModel, probabilityMatrix, approxRowGroupIndices, rowSubstitutions, stateActionRewards, rewardEntryToEvalTableMapping, &constantRewEntry);
rewardModels.insert(std::pair<std::string, storm::models::sparse::StandardRewardModel<ConstantType>>("", storm::models::sparse::StandardRewardModel<ConstantType>(boost::optional<std::vector<ConstantType>>(), boost::optional<std::vector<ConstantType>>(std::move(stateActionRewards)))));
}
//Obtain other model ingredients and the approximation model itself
storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;
switch(parametricModel.getType()){
case storm::models::ModelType::Dtmc:
this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
break;
case storm::models::ModelType::Mdp:
//TODO: use a two-player-game here
this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
break;
default:
STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Tried to build an Approximation model for an unsupported model type");
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<std::size_t> newIndices(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates()); //initialize with some illegal index
std::size_t newIndex=0;
for(auto const& index : maybeStates){
newIndices[index]=newIndex;
++newIndex;
}
//translate the matrixEntryToEvalTableMapping into the actual probability mapping
typename storm::storage::SparseMatrix<ConstantType>::index_type matrixRow=0;
auto tableEntry=matrixEntryToEvalTableMapping.begin();
for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
for (; matrixRow<approxRowGroupIndices[row+1]; ++matrixRow){
auto approxModelEntry = this->model->getTransitionMatrix().getRow(matrixRow).begin();
for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){
if(*tableEntry == &constantProbEntry){
STORM_LOG_THROW(storm::utility::isConstant(parEntry.getValue()), storm::exceptions::UnexpectedException, "Expected a constant matrix entry but got a non-constant one.");
approxModelEntry->setValue(storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parEntry.getValue())));
} else {
this->probabilityMapping.emplace_back(std::make_pair(&((*tableEntry)->second), approxModelEntry));
}
++approxModelEntry;
++tableEntry;
}
}
}
//Now pre-compute the information for the equation system.
std::vector<std::size_t> rowSubstitutions;// the substitution used in every row
initializeProbabilities(parametricModel, newIndices, rowSubstitutions);
if(this->computeRewards){
//the same for rewards
auto approxModelRewardEntry = this->model->getUniqueRewardModel()->second.getStateActionRewardVector().begin();
for (auto tableEntry : rewardEntryToEvalTableMapping){
if(tableEntry != &constantRewEntry){
//insert pointer to minValue, pointer to maxValue, iterator to reward vector entry
this->rewardMapping.emplace_back(std::make_tuple(&(tableEntry->second.first), &(tableEntry->second.second), approxModelRewardEntry));
}
++approxModelRewardEntry;
}
//Get the sets of reward parameters that map to "CHOSEOPTIMAL".
this->choseOptimalRewardPars = std::vector<std::set<VariableType>>(this->rewardSubstitutions.size());
for(std::size_t index = 0; index<this->rewardSubstitutions.size(); ++index){
for(auto const& sub : this->rewardSubstitutions[index]){
if(sub.second==TypeOfBound::CHOSEOPTIMAL){
this->choseOptimalRewardPars[index].insert(sub.first);
}
}
initializeRewards(parametricModel, newIndices, rowSubstitutions);
}
this->matrixData.assignment.shrink_to_fit();
this->vectorData.assignment.shrink_to_fit();
if(parametricModel.getType()==storm::models::ModelType::Mdp){
initializePlayer1Matrix(parametricModel, newIndices);
}
this->eqSysResult = std::vector<ConstantType>(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one<ConstantType>() : ConstantType(0.5));
this->eqSysInitIndex = newIndices[initialState];
}
template<typename ParametricSparseModelType, typename ConstantType>
void ApproximationModel<ParametricSparseModelType, ConstantType>::initializeProbabilities(ParametricSparseModelType const& parametricModel,
storm::storage::SparseMatrix<ConstantType>& probabilityMatrix,
std::vector<typename storm::storage::SparseMatrix<ConstantType>::index_type>& approxRowGroupIndices,
std::vector<std::size_t>& rowSubstitutions,
std::vector<ProbTableEntry*>& matrixEntryToEvalTableMapping,
ProbTableEntry* constantEntry) {
void ApproximationModel<ParametricSparseModelType, ConstantType>::initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t>& rowSubstitutions) {
STORM_LOG_DEBUG("Approximation model initialization for probabilities");
// Run through the rows of the original model and obtain the different substitutions, the probability evaluation table,
// a probability matrix with some dummy entries, and the mapping between the two
//TODO: if the parametricModel is nondeterministic, we will need a matrix for a two player game. For now, we assume the two players cooperate, i.e. we get a simple MDP
bool customRowGroups = parametricModel.getTransitionMatrix().hasNontrivialRowGrouping();
auto curParamRowGroup = parametricModel.getTransitionMatrix().getRowGroupIndices().begin();
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(0, //unknown number of rows
parametricModel.getTransitionMatrix().getColumnCount(),
/* First run: get a matrix with dummy entries at the new positions.
* This matrix will have a row group for every row in the original matrix,
* each rowgroup containing 2^#par rows, where #par is the number of parameters that occur in the original row.
* We also store the substitution that needs to be applied for each row.
*/
ConstantType dummyValue = storm::utility::one<ConstantType>();
auto numOfMaybeStates = this->maybeStates.getNumberOfSetBits();
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(numOfMaybeStates, //exact number of rows is unknown at this point, but at least this many
numOfMaybeStates, //columns
0, //Unknown number of entries
true, //force dimensions
false, // no force dimensions
true, //will have custom row grouping
parametricModel.getNumberOfStates());
approxRowGroupIndices.reserve(parametricModel.getTransitionMatrix().getRowCount()+1);
std::size_t numOfNonConstEntries=0;
typename storm::storage::SparseMatrix<ConstantType>::index_type apprModelMatrixRow=0;
for(typename storm::storage::SparseMatrix<ParametricType>::index_type paramRow=0; paramRow < parametricModel.getTransitionMatrix().getRowCount(); ++paramRow ){
approxRowGroupIndices.push_back(apprModelMatrixRow);
if(customRowGroups && paramRow==*curParamRowGroup){
//Only make a new row group if the parametric model matrix has a new row group at this position
matrixBuilder.newRowGroup(apprModelMatrixRow);
++curParamRowGroup;
} else if (!customRowGroups){
matrixBuilder.newRowGroup(apprModelMatrixRow);
}
numOfMaybeStates); //exact number of rowgroups is unknown at this point, but at least this many
rowSubstitutions.push_back(numOfMaybeStates);
std::size_t curRow = 0;
for (auto oldRowGroup : this->maybeStates){
for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){
matrixBuilder.newRowGroup(curRow);
// Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper}
std::set<VariableType> occurringVariables;
for(auto const& entry : parametricModel.getTransitionMatrix().getRow(paramRow)){
storm::utility::region::gatherOccurringVariables(entry.getValue(), occurringVariables);
for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){
if(this->maybeStates.get(oldEntry.getColumn())){
storm::utility::region::gatherOccurringVariables(oldEntry.getValue(), occurringVariables);
}
}
std::size_t numOfSubstitutions=1ull<<occurringVariables.size(); //=2^(#variables). Note that there is still 1 substitution when #variables==0 (i.e.,the empty substitution)
std::size_t numOfSubstitutions=1ull<<occurringVariables.size(); //=2^(#variables). Note that there is still 1 substitution when #variables==0 (the empty substitution)
for(uint_fast64_t substitutionId=0; substitutionId<numOfSubstitutions; ++substitutionId){
//compute actual substitution from substitutionId by interpreting the Id as a bit sequence.
//the occurringVariables.size() least significant bits of substitutionId will always represent the substitution that we have to consider
//the occurringVariables.size() least significant bits of substitutionId will represent the substitution that we have to consider
//(00...0 = lower bounds for all parameters, 11...1 = upper bounds for all parameters)
std::map<VariableType, TypeOfBound> currSubstitution;
std::size_t parameterIndex=0;
std::size_t parameterIndex=0ull;
for(auto const& parameter : occurringVariables){
if((substitutionId>>parameterIndex)%2==0){
currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER));
@ -162,92 +114,129 @@ namespace storm {
}
++parameterIndex;
}
std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->probabilitySubstitutions, std::move(currSubstitution));
std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->funcSubData.substitutions, std::move(currSubstitution));
rowSubstitutions.push_back(substitutionIndex);
//For every substitution, run again through the row and add an entry in matrixEntryToEvalTableMapping as well as dummy entries in the matrix
//For every substitution, run again through the row and add a dummy entry
//Note that this is still executed once, even if no parameters occur.
ConstantType dummyEntry=storm::utility::one<ConstantType>();
for(auto const& entry : parametricModel.getTransitionMatrix().getRow(paramRow)){
if(storm::utility::isConstant(entry.getValue())){
matrixEntryToEvalTableMapping.emplace_back(constantEntry);
for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){
if(this->maybeStates.get(oldEntry.getColumn())){
matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyValue);
}
}
}
++curRow;
}
}
this->matrixData.matrix=matrixBuilder.build();
//Now run again through both matrices to get the remaining ingredients of the matrixData and vectorData
this->matrixData.assignment.reserve(this->matrixData.matrix.getEntryCount());
this->vectorData.vector = std::vector<ConstantType>(this->matrixData.matrix.getRowCount()); //Important to initialize here since iterators have to remain valid
auto vectorIt = this->vectorData.vector.begin();
this->vectorData.assignment.reserve(vectorData.vector.size());
std::size_t curRowGroup = 0;
for(auto oldRowGroup : this->maybeStates){
for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){
ParametricType targetProbability = storm::utility::region::getNewFunction<ParametricType, CoefficientType>(storm::utility::zero<CoefficientType>());
if(!this->computeRewards){
for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){
if(this->targetStates.get(oldEntry.getColumn())){
targetProbability += oldEntry.getValue();
}
}
}
//Recall: Every row in the old matrix has a row group in the newly created one.
//We will now run through every row that belongs to the rowGroup associated with oldRow.
for (curRow = this->matrixData.matrix.getRowGroupIndices()[curRowGroup]; curRow < this->matrixData.matrix.getRowGroupIndices()[curRowGroup+1]; ++curRow){
auto eqSysMatrixEntry = this->matrixData.matrix.getRow(curRow).begin();
for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){
if(this->maybeStates.get(oldEntry.getColumn())){
STORM_LOG_THROW(eqSysMatrixEntry->getColumn()==newIndices[oldEntry.getColumn()], storm::exceptions::UnexpectedException, "old and new entries do not match");
if(storm::utility::isConstant(oldEntry.getValue())){
eqSysMatrixEntry->setValue(storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(oldEntry.getValue())));
} else {
auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(oldEntry.getValue(), rowSubstitutions[curRow]), 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;
}
}
if(!this->computeRewards){
if(storm::utility::isConstant(targetProbability)){
*vectorIt = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(targetProbability));
} else {
++numOfNonConstEntries;
auto evalTableIt = this->probabilityEvaluationTable.insert(ProbTableEntry(FunctionSubstitution(entry.getValue(), substitutionIndex), storm::utility::zero<ConstantType>())).first;
matrixEntryToEvalTableMapping.emplace_back(&(*evalTableIt));
auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(targetProbability, rowSubstitutions[curRow]), dummyValue)).first;
this->vectorData.assignment.emplace_back(std::make_pair(vectorIt, &(functionsIt->second)));
*vectorIt = dummyValue;
}
matrixBuilder.addNextValue(apprModelMatrixRow, entry.getColumn(), dummyEntry);
dummyEntry=storm::utility::zero<ConstantType>();
}
++apprModelMatrixRow;
++vectorIt;
}
++curRowGroup;
}
this->probabilityMapping.reserve(numOfNonConstEntries);
probabilityMatrix=matrixBuilder.build();
approxRowGroupIndices.push_back(apprModelMatrixRow);
}
STORM_LOG_THROW(vectorIt==this->vectorData.vector.end(), storm::exceptions::UnexpectedException, "initProbs: The size of the eq-sys vector is not as it was expected");
this->matrixData.matrix.updateNonzeroEntryCount();
}
template<typename ParametricSparseModelType, typename ConstantType>
void ApproximationModel<ParametricSparseModelType, ConstantType>::initializeRewards(ParametricSparseModelType const& parametricModel,
storm::storage::SparseMatrix<ConstantType> const& probabilityMatrix,
std::vector<typename storm::storage::SparseMatrix<ConstantType>::index_type> const& approxRowGroupIndices,
std::vector<std::size_t> const& rowSubstitutions,
std::vector<ConstantType>& stateActionRewardVector,
std::vector<RewTableEntry*>& rewardEntryToEvalTableMapping,
RewTableEntry* constantEntry){
void ApproximationModel<ParametricSparseModelType, ConstantType>::initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t> const& rowSubstitutions){
STORM_LOG_DEBUG("Approximation model initialization for Rewards");
STORM_LOG_THROW(parametricModel.getType()==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Rewards are only supported for DTMCs (yet)"); // Todo : check if this code also works for rewards on mdps (we should then consider state action rewards of the original model and approxRowGroupIndices)
STORM_LOG_THROW(parametricModel.getType()==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Rewards are only supported for DTMCs (yet)");
//Note: Since the original model is assumed to be a DTMC, there is no outgoing transition of a maybeState that leads to an infinity state.
//Hence, we do not have to set entries of the eqSys vector to infinity (as it would be required for mdp model checking...)
STORM_LOG_THROW(this->vectorData.vector.size()==this->matrixData.matrix.getRowCount(), storm::exceptions::UnexpectedException, "The size of the eq-sys vector does not match to the number of rows in the eq-sys matrix");
this->vectorData.assignment.reserve(vectorData.vector.size());
// run through the state reward vector of the parametric model.
// Constant entries can be set directly.
// For Parametric entries we set a dummy value and insert one entry to the rewardEntryEvalTableMapping
stateActionRewardVector.reserve(probabilityMatrix.getRowCount());
rewardEntryToEvalTableMapping.reserve(probabilityMatrix.getRowCount());
std::size_t numOfNonConstRewEntries=0;
for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
std::set<VariableType> occurringRewVariables;
std::set<VariableType> occurringProbVariables;
if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])){
ConstantType reward = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state]));
//Add one of these entries for every row in the row group of state
for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow<probabilityMatrix.getRowGroupIndices()[state+1]; ++matrixRow){
stateActionRewardVector.emplace_back(reward);
rewardEntryToEvalTableMapping.emplace_back(constantEntry);
// For Parametric entries we set a dummy value and insert the corresponding function and the assignment
ConstantType dummyValue = storm::utility::one<ConstantType>();
auto vectorIt = this->vectorData.vector.begin();
for(auto oldState : this->maybeStates){
if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState])){
ConstantType reward = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState]));
//Add one of these entries for every row in the row group of oldState
for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[state+1]; ++matrixRow){
*vectorIt = reward;
++vectorIt;
}
} else {
storm::utility::region::gatherOccurringVariables(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], occurringRewVariables);
//For each row in the row group of state, we get the corresponding substitution and tableIndex
// We might find out that the reward is independent of the probability variables (and will thus be independent of nondeterministic choices)
// In that case, the reward function and the substitution will not change and thus we can use the same table index
std::set<VariableType> occurringRewVariables;
storm::utility::region::gatherOccurringVariables(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState], occurringRewVariables);
// For each row in the row group of oldState, we get the corresponding substitution and insert the FunctionSubstitution
// We might find out that the reward is independent of the probability parameters (and will thus be independent of nondeterministic choices)
// In that case, the reward function and the substitution will not change and thus we can use the same FunctionSubstitution
bool rewardDependsOnProbVars=true;
RewTableEntry* evalTablePtr;
for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow<probabilityMatrix.getRowGroupIndices()[state+1]; ++matrixRow){
std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::iterator functionsIt;
for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[oldState+1]; ++matrixRow){
auto probabilitySub = this->funcSubData.substitutions[rowSubstitutions[matrixRow]];
if(rewardDependsOnProbVars){ //always executed in first iteration
rewardDependsOnProbVars=false; //Assume that independent...
std::map<VariableType, TypeOfBound> rewardSubstitution;
//Get the correct substitution for this matrixRow
std::map<VariableType, TypeOfBound> substitution;
for(auto const& rewardVar : occurringRewVariables){
typename std::map<VariableType, TypeOfBound>::iterator const substitutionIt=this->probabilitySubstitutions[rowSubstitutions[matrixRow]].find(rewardVar);
if(substitutionIt==this->probabilitySubstitutions[rowSubstitutions[matrixRow]].end()){
rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
typename std::map<VariableType, TypeOfBound>::iterator const substitutionIt = probabilitySub.find(rewardVar);
if(substitutionIt==probabilitySub.end()){
substitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
} else {
rewardSubstitution.insert(*substitutionIt);
substitution.insert(*substitutionIt);
rewardDependsOnProbVars=true; //.. assumption wrong
}
}
std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->rewardSubstitutions, std::move(rewardSubstitution));
auto evalTableIt = this->rewardEvaluationTable.insert(
RewTableEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], substitutionIndex),
std::pair<ConstantType, ConstantType>(storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>()))
).first;
evalTablePtr=&(*evalTableIt);
// insert the substitution and the FunctionSubstitution
std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->funcSubData.substitutions, std::move(substitution));
functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], substitutionIndex), dummyValue)).first;
}
//insert table entries and dummy data
stateActionRewardVector.emplace_back(storm::utility::zero<ConstantType>());
rewardEntryToEvalTableMapping.emplace_back(evalTablePtr);
++numOfNonConstRewEntries;
//insert assignment and dummy data
this->vectorData.assignment.emplace_back(std::make_pair(vectorIt, &(functionsIt->second)));
*vectorIt = dummyValue;
++vectorIt;
}
}
}
this->rewardMapping.reserve(numOfNonConstRewEntries);
STORM_LOG_THROW(vectorIt==this->vectorData.vector.end(), storm::exceptions::UnexpectedException, "initRewards: The size of the eq-sys vector is not as it was expected");
}
@ -258,136 +247,100 @@ namespace storm {
}
template<typename ParametricSparseModelType, typename ConstantType>
std::shared_ptr<storm::models::sparse::Model<ConstantType>> const& ApproximationModel<ParametricSparseModelType, ConstantType>::getModel() const {
return this->model;
std::vector<ConstantType> ApproximationModel<ParametricSparseModelType, ConstantType>::computeValues(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
instantiate(region, computeLowerBounds);
invokeSolver(computeLowerBounds);
std::vector<ConstantType> 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<ConstantType>() : storm::utility::one<ConstantType>());
storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>());
return result;
}
template<typename ParametricSparseModelType, typename ConstantType>
void ApproximationModel<ParametricSparseModelType, ConstantType>::instantiate(const ParameterRegion<ParametricType>& region) {
//Instantiate the substitutions
std::vector<std::map<VariableType, CoefficientType>> instantiatedSubs(this->probabilitySubstitutions.size());
for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->probabilitySubstitutions.size(); ++substitutionIndex){
for(std::pair<VariableType, TypeOfBound> const& sub : this->probabilitySubstitutions[substitutionIndex]){
switch(sub.second){
case TypeOfBound::LOWER:
instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first)));
break;
case TypeOfBound::UPPER:
instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getUpperBound(sub.first)));
break;
default:
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound");
}
}
}
//write entries into evaluation table
for(auto& tableEntry : this->probabilityEvaluationTable){
tableEntry.second=storm::utility::region::convertNumber<ConstantType>(
storm::utility::region::evaluateFunction(
tableEntry.first.getFunction(),
instantiatedSubs[tableEntry.first.getSubstitution()]
)
);
}
//write the instantiated values to the matrix according to the mapping
for(auto& mappingPair : this->probabilityMapping){
mappingPair.second->setValue(*(mappingPair.first));
ConstantType ApproximationModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
instantiate(region, computeLowerBounds);
invokeSolver(computeLowerBounds);
return this->eqSysResult[this->eqSysInitIndex];
}
if(this->computeRewards){
template<typename ParametricSparseModelType, typename ConstantType>
void ApproximationModel<ParametricSparseModelType, ConstantType>::instantiate(const ParameterRegion<ParametricType>& region, bool computeLowerBounds) {
//Instantiate the substitutions
std::vector<std::map<VariableType, CoefficientType>> instantiatedRewardSubs(this->rewardSubstitutions.size());
for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->rewardSubstitutions.size(); ++substitutionIndex){
for(std::pair<VariableType, TypeOfBound> const& sub : this->rewardSubstitutions[substitutionIndex]){
std::vector<std::map<VariableType, CoefficientType>> instantiatedSubs(this->funcSubData.substitutions.size());
std::vector<std::set<VariableType>> choseOptimalParameters(this->funcSubData.substitutions.size());
for(std::size_t substitutionIndex=0; substitutionIndex<this->funcSubData.substitutions.size(); ++substitutionIndex){
for(std::pair<VariableType, TypeOfBound> const& sub : this->funcSubData.substitutions[substitutionIndex]){
switch(sub.second){
case TypeOfBound::LOWER:
instantiatedRewardSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first)));
instantiatedSubs[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)));
instantiatedSubs[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<ConstantType>()));
instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, storm::utility::one<CoefficientType>()));
choseOptimalParameters[substitutionIndex].insert(sub.first);
break;
default:
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound when instantiating a reward substitution. Index: " << substitutionIndex << " TypeOfBound: "<< ((int)sub.second));
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound");
}
}
}
//write entries into evaluation table
for(auto& tableEntry : this->rewardEvaluationTable){
auto& funcSub = tableEntry.first;
auto& minMax = tableEntry.second;
minMax.first=storm::utility::infinity<ConstantType>();
minMax.second=storm::utility::zero<ConstantType>();
//write function+substitution results into placeholders
for(auto& functionResult : this->funcSubData.functions){
auto& funcSub = functionResult.first;
auto& result = functionResult.second;
result = computeLowerBounds ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>();
//Iterate over the different combinations of lower and upper bounds and update the min and max values
auto const& vertices=region.getVerticesOfRegion(this->choseOptimalRewardPars[funcSub.getSubstitution()]);
auto const& vertices=region.getVerticesOfRegion(this->choseOptimalParameters[funcSub.getSubstitution()]);
for(auto const& vertex : vertices){
//extend the substitution
for(auto const& vertexSub : vertex){
instantiatedRewardSubs[funcSub.getSubstitution()][vertexSub.first]=vertexSub.second;
instantiatedSubs[funcSub.getSubstitution()][vertexSub.first]=vertexSub.second;
}
//evaluate the function
ConstantType currValue = storm::utility::region::convertNumber<ConstantType>(
storm::utility::region::evaluateFunction(
funcSub.getFunction(),
instantiatedRewardSubs[funcSub.getSubstitution()]
instantiatedSubs[funcSub.getSubstitution()]
)
);
minMax.first=std::min(minMax.first, currValue);
minMax.second=std::max(minMax.second, currValue);
}
}
//Note: the rewards are written to the model as soon as the optimality type is known (i.e. in computeValues)
result = computeLowerBounds ? std::min(result, currValue) : std::max(result, currValue);
}
}
template<typename ParametricSparseModelType, typename ConstantType>
std::unique_ptr<storm::modelchecker::CheckResult> ApproximationModel<ParametricSparseModelType, ConstantType>::computeValues(storm::solver::OptimizationDirection const& approximationOpDir) {
std::unique_ptr<storm::modelchecker::AbstractModelChecker> modelChecker;
switch(this->getModel()->getType()){
case storm::models::ModelType::Mdp:
modelChecker = std::make_unique<storm::modelchecker::SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<ConstantType>>>(*this->model->template as<storm::models::sparse::Mdp<ConstantType>>());
break;
case storm::models::ModelType::S2pg:
//TODO: instantiate right model checker here
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Approximation with two player games not yet implemented");
break;
default:
STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Tried to build an approximation model for an unsupported model type");
//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));
}
std::unique_ptr<storm::modelchecker::CheckResult> resultPtr;
//TODO: use this when two player games are available
boost::optional<storm::solver::OptimizationDirection> formulaOpDir = storm::logic::isLowerBound(this->formula->getComparisonType()) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize;
if(this->computeRewards){
//write the reward values into the model
switch(approximationOpDir){
case storm::solver::OptimizationDirection::Minimize:
for(auto& mappingTriple : this->rewardMapping){
*std::get<2>(mappingTriple)=*std::get<0>(mappingTriple);
for(auto& assignment : this->vectorData.assignment){
*assignment.first=*assignment.second;
}
break;
case storm::solver::OptimizationDirection::Maximize:
for(auto& mappingTriple : this->rewardMapping){
*std::get<2>(mappingTriple)=*std::get<1>(mappingTriple);
}
break;
default:
STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given optimality Type is not supported.");
}
//perform model checking
boost::optional<std::string> noRewardModelName; //it should be uniquely given
resultPtr = modelChecker->computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula(), noRewardModelName, false, approximationOpDir);
}
else {
//perform model checking
resultPtr = modelChecker->computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula(), false, approximationOpDir);
template<>
void ApproximationModel<storm::models::sparse::Dtmc<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds){
storm::solver::SolveGoal goal(computeLowerBounds);
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> solver = storm::solver::configureMinMaxLinearEquationSolver(goal, storm::utility::solver::MinMaxLinearEquationSolverFactory<double>(), this->matrixData.matrix);
solver->solveEquationSystem(this->eqSysResult, this->vectorData.vector);
}
return resultPtr;
template<>
void ApproximationModel<storm::models::sparse::Mdp<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds){
storm::solver::SolveGoal player2Goal(computeLowerBounds);
std::unique_ptr<storm::solver::GameSolver<double>> solver = storm::utility::solver::GameSolverFactory<double>().create(this->player1Matrix, this->matrixData.matrix);
solver->solveGame(this->player1Goal.direction(), player2Goal.direction(), this->eqSysResult, this->vectorData.vector);
}
#ifdef STORM_HAVE_CARL
template class ApproximationModel<storm::models::sparse::Dtmc<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
template class ApproximationModel<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;

106
src/modelchecker/region/ApproximationModel.h

@ -17,6 +17,7 @@
#include "src/logic/Formulas.h"
#include "src/models/sparse/Model.h"
#include "src/storage/SparseMatrix.h"
#include "src/solver/SolveGoal.h"
namespace storm {
namespace modelchecker {
@ -30,28 +31,28 @@ namespace storm {
/*!
* Creates an Approximation 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)
* @note This will not check whether approximation is applicable
*/
ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula);
virtual ~ApproximationModel();
/*!
* returns the underlying model
* Instantiates the approximation model w.r.t. the given region.
* Then computes and returns the approximated reachability probabilities or reward values for every state.
* If computeLowerBounds is true, the computed values will be a lower bound for the actual values. Otherwise, we get upper bounds,
*/
std::shared_ptr<storm::models::sparse::Model<ConstantType>> const& getModel() const;
std::vector<ConstantType> computeValues(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
/*!
* Instantiates the underlying model according to the given region
* Instantiates the approximation model w.r.t. the given region.
* Then computes and returns the approximated reachability probabilities or reward value for the initial state.
* If computeLowerBounds is true, the computed value will be a lower bound for the actual value. Otherwise, we get an upper bound.
*/
void instantiate(ParameterRegion<ParametricType> const& region);
/*!
* Returns the approximated reachability probabilities or reward values for every state.
* Undefined behavior if model has not been instantiated first!
* @param approximationOpDir Use MAXIMIZE to get upper bounds or MINIMIZE to get lower bounds
*/
std::unique_ptr<storm::modelchecker::CheckResult> computeValues(storm::solver::OptimizationDirection const& approximationOpDir);
ConstantType computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
private:
//This enum helps to store how a parameter will be substituted.
@ -121,42 +122,59 @@ namespace storm {
}
};
typedef typename std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::value_type ProbTableEntry;
typedef typename std::unordered_map<FunctionSubstitution, std::pair<ConstantType, ConstantType>, FuncSubHash>::value_type RewTableEntry;
void initializeProbabilities(ParametricSparseModelType const& parametricModel, storm::storage::SparseMatrix<ConstantType>& probabilityMatrix, std::vector<typename storm::storage::SparseMatrix<ConstantType>::index_type>& approxRowGroupIndices, std::vector<std::size_t>& rowSubstitutions, std::vector<ProbTableEntry*>& matrixEntryToEvalTableMapping, ProbTableEntry* constantEntry);
void initializeRewards(ParametricSparseModelType const& parametricModel, storm::storage::SparseMatrix<ConstantType> const& probabilityMatrix, std::vector<typename storm::storage::SparseMatrix<ConstantType>::index_type> const& approxRowGroupIndices, std::vector<std::size_t> const& rowSubstitutions, std::vector<ConstantType>& stateActionRewardVector, std::vector<RewTableEntry*>& rewardEntryToEvalTableMapping, RewTableEntry* constantEntry);
//The Model with which we work
std::shared_ptr<storm::models::sparse::Model<ConstantType>> model;
//The formula for which we will compute the values
std::shared_ptr<storm::logic::OperatorFormula> formula;
typedef typename std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::value_type FunctionEntry;
// typedef typename std::unordered_map<FunctionSubstitution, std::pair<ConstantType, ConstantType>, FuncSubHash>::value_type RewTableEntry;
void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t>& rowSubstitutions);
void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t> const& rowSubstitutions);
void initializePlayer1Matrix(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
void instantiate(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
void invokeSolver(bool computeLowerBounds);
//Some designated states in the original 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<ConstantType> 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;
//Player 1 represents the nondeterminism of the given mdp (so, this is irrelevant if we approximate values of a DTMC)
storm::solver::SolveGoal player1Goal;
storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix;
// We store one (unique) entry for every occurring pair of a non-constant function and substitution.
// Whenever a region is given, we can then evaluate the functions (w.r.t. the corresponding substitutions)
// and store the result to the target value of this map
std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash> probabilityEvaluationTable;
//For rewards, we map to the minimal value and the maximal value (depending on the CHOSEOPTIMAL parameters).
std::unordered_map<FunctionSubstitution, std::pair<ConstantType, ConstantType>, FuncSubHash> rewardEvaluationTable;
/* 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 pair of a non-constant function and substitution.
* The map points to some ConstantType value which serves as placeholder.
* When instantiating the model, the evaluated result of every function + substitution is stored in the corresponding placeholder.
* For rewards, however, we might need a minimal and a maximal value which is why there are two paceholders.
* There is an assignment that connects every non-constant matrix (or: vector) 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 FuncSubData{
// the occurring probability functions together with the corresponding placeholders for the result
std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash> functions;
// the occurring reward functions together with the corresponding placeholders for the minimal and maximal result
// //std::unordered_map<FunctionSubstitution, std::pair<ConstantType, ConstantType>, FuncSubHash> rewFunctions;
//Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region)
std::vector<std::map<VariableType, TypeOfBound>> 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 (i.e. CHOSEOPTIMAL parameters)
std::vector<std::map<VariableType, TypeOfBound>> rewardSubstitutions;
std::vector<std::set<VariableType>> choseOptimalRewardPars;
//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<std::pair<ConstantType*, typename storm::storage::SparseMatrix<ConstantType>::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.
std::vector<std::tuple<ConstantType*, ConstantType*, typename std::vector<ConstantType>::iterator>> rewardMapping;
std::vector<std::map<VariableType, TypeOfBound>> substitutions;
//For rewards, we also store the parameters for which the correct substitution depends on the region and whether to minimize/maximize (i.e. TypeOfBound==CHOSEOPTIMAL)
// // std::vector<std::map<VariableType, TypeOfBound>> rewSubs;
// //std::vector<std::set<VariableType>> choseOptimalRewardPars;
} funcSubData;
struct MatrixData {
storm::storage::SparseMatrix<ConstantType> matrix; //The matrix itself.
std::vector<std::pair<typename storm::storage::SparseMatrix<ConstantType>::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders
// std::vector<std::tuple<ConstantType*, ConstantType*, typename std::vector<ConstantType>::iterator>> rewardMapping;
} matrixData;
struct VectorData {
std::vector<ConstantType> vector; //The vector itself.
std::vector<std::pair<typename std::vector<ConstantType>::iterator, ConstantType*>> assignment; // Connection of vector entries with placeholders
} vectorData;
};
} //namespace region

106
src/modelchecker/region/SamplingModel.cpp

@ -10,12 +10,12 @@
#include "src/models/sparse/Dtmc.h"
#include "src/models/sparse/Mdp.h"
#include "src/models/ModelType.h"
#include "models/sparse/StandardRewardModel.h"
#include "src/utility/solver.h"
#include "src/models/sparse/StandardRewardModel.h"
#include "src/solver/LinearEquationSolver.h"
#include "src/solver/MinMaxLinearEquationSolver.h"
#include "src/utility/macros.h"
#include "src/utility/region.h"
#include "src/utility/solver.h"
#include "src/utility/vector.h"
#include "src/exceptions/UnexpectedException.h"
#include "src/exceptions/InvalidArgumentException.h"
@ -65,24 +65,21 @@ namespace storm {
this->eqSysResult = std::vector<ConstantType>(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one<ConstantType>() : ConstantType(0.5));
this->eqSysInitIndex = newIndices[initialState];
this->solveGoal = storm::solver::SolveGoal(storm::logic::isLowerBound(formula->getComparisonType()));
}
template<typename ParametricSparseModelType, typename ConstantType>
void SamplingModel<ParametricSparseModelType, ConstantType>::initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices){
//First run: get a matrix with dummy entries at the new positions and fill the vectorData
//First run: get a matrix with dummy entries at the new positions
ConstantType dummyValue = storm::utility::one<ConstantType>();
bool addRowGroups = parametricModel.getTransitionMatrix().hasNontrivialRowGrouping();
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();
auto numOfMaybeStates = this->maybeStates.getNumberOfSetBits();
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(addRowGroups ? parametricModel.getTransitionMatrix().getRowCount() : numOfMaybeStates,
numOfMaybeStates,
parametricModel.getTransitionMatrix().getEntryCount(),
false, // no force dimensions
addRowGroups,
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){
@ -90,7 +87,6 @@ namespace storm {
}
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<ParametricType, CoefficientType>(storm::utility::zero<CoefficientType>());
for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){
if(this->maybeStates.get(oldEntry.getColumn())){
//check if we need to insert a diagonal entry
@ -105,36 +101,28 @@ namespace storm {
//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<ConstantType>());
}
if(!this->computeRewards){
if(storm::utility::isConstant(targetProbability)){
this->vectorData.vector.push_back(storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(targetProbability)));
} else {
typename std::unordered_map<ParametricType, ConstantType>::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;
}
}
this->matrixData.matrix=matrixBuilder.build();
//Now run again through both matrices to get the remaining ingredients of the matrixData.
//Now run again through both matrices to get the remaining ingredients of the matrixData and vectorData.
//Note that we need the matrix (I-P) in case of a dtmc.
this->matrixData.assignment.reserve(this->matrixData.matrix.getEntryCount());
this->vectorData.vector = std::vector<ConstantType>(this->matrixData.matrix.getRowCount()); //Important to initialize here since iterators have to remain valid
auto vectorIt = this->vectorData.vector.begin();
this->vectorData.assignment.reserve(vectorData.vector.size());
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();
ParametricType targetProbability = storm::utility::region::getNewFunction<ParametricType, CoefficientType>(storm::utility::zero<CoefficientType>());
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()]){
@ -157,40 +145,58 @@ namespace storm {
typename std::unordered_map<ParametricType, ConstantType>::iterator functionsIt;
if(isDtmc){
if(eqSysMatrixEntry->getColumn()==newIndices[oldRowGroup]){ //Diagonal entries get 1-f(x)
functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::one<ParametricType>()-oldEntry.getValue(), dummyValue)).first;
functionsIt = this->functions.insert(FunctionEntry(storm::utility::one<ParametricType>()-oldEntry.getValue(), dummyValue)).first;
} else {
functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::zero<ParametricType>()-oldEntry.getValue(), dummyValue)).first;
functionsIt = this->functions.insert(FunctionEntry(storm::utility::zero<ParametricType>()-oldEntry.getValue(), dummyValue)).first;
}
} else {
functionsIt = this->matrixData.functions.insert(FunctionEntry(oldEntry.getValue(), dummyValue)).first;
functionsIt = this->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;
}
else if(!this->computeRewards && this->targetStates.get(oldEntry.getColumn())){
targetProbability += oldEntry.getValue();
}
}
if(!this->computeRewards){
if(storm::utility::isConstant(targetProbability)){
*vectorIt = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(targetProbability));
} else {
typename std::unordered_map<ParametricType, ConstantType>::iterator functionsIt = this->functions.insert(FunctionEntry(targetProbability, dummyValue)).first;
this->vectorData.assignment.emplace_back(std::make_pair(vectorIt, &(functionsIt->second)));
*vectorIt = dummyValue;
}
}
++vectorIt;
++curRow;
}
}
STORM_LOG_THROW(vectorIt==this->vectorData.vector.end(), storm::exceptions::UnexpectedException, "initProbs: The size of the eq-sys vector is not as it was expected");
this->matrixData.matrix.updateNonzeroEntryCount();
}
template<typename ParametricSparseModelType, typename ConstantType>
void SamplingModel<ParametricSparseModelType, ConstantType>::initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> 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());
//Run through the state reward vector... Note that this only works for dtmcs
STORM_LOG_THROW(parametricModel.getType()==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Rewards are only supported for DTMCs (yet)");
STORM_LOG_THROW(this->vectorData.vector.size()==this->matrixData.matrix.getRowCount(), storm::exceptions::UnexpectedException, "The size of the eq-sys vector does not match to the number of rows in the eq-sys matrix");
this->vectorData.assignment.reserve(vectorData.vector.size());
ConstantType dummyValue = storm::utility::one<ConstantType>();
auto vectorIt = this->vectorData.vector.begin();
for(auto state : this->maybeStates){
if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])){
this->vectorData.vector.push_back(storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])));
*vectorIt = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state]));
} else {
typename std::unordered_map<ParametricType, ConstantType>::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);
typename std::unordered_map<ParametricType, ConstantType>::iterator functionsIt = this->functions.insert(FunctionEntry(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], dummyValue)).first;
this->vectorData.assignment.emplace_back(std::make_pair(vectorIt, &(functionsIt->second)));
*vectorIt = dummyValue;
}
++vectorIt;
}
STORM_LOG_THROW(vectorIt==this->vectorData.vector.end(), storm::exceptions::UnexpectedException, "initRewards: The size of the eq-sys vector is not as it was expected");
}
template<typename ParametricSparseModelType, typename ConstantType>
@ -198,14 +204,31 @@ namespace storm {
//Intentionally left empty
}
template<typename ParametricSparseModelType, typename ConstantType>
std::vector<ConstantType> SamplingModel<ParametricSparseModelType, ConstantType>::computeValues() {
invokeSolver();
std::vector<ConstantType> 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<ConstantType>() : storm::utility::one<ConstantType>());
storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>());
return result;
}
template<typename ParametricSparseModelType, typename ConstantType>
ConstantType SamplingModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue() {
invokeSolver();
return this->eqSysResult[this->eqSysInitIndex];
}
template<typename ParametricSparseModelType, typename ConstantType>
void SamplingModel<ParametricSparseModelType, ConstantType>::instantiate(std::map<VariableType, CoefficientType>const& point) {
//write results into the placeholders
for(auto& functionResult : this->matrixData.functions){
for(auto& functionResult : this->functions){
functionResult.second=storm::utility::region::convertNumber<ConstantType>(
storm::utility::region::evaluateFunction(functionResult.first, point));
}
for(auto& functionResult : this->vectorData.functions){
for(auto& functionResult : this->functions){
functionResult.second=storm::utility::region::convertNumber<ConstantType>(
storm::utility::region::evaluateFunction(functionResult.first, point));
}
@ -219,23 +242,6 @@ namespace storm {
}
}
template<typename ParametricSparseModelType, typename ConstantType>
std::vector<ConstantType> SamplingModel<ParametricSparseModelType, ConstantType>::computeValues() {
invokeSolver();
std::vector<ConstantType> 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<ConstantType>() : storm::utility::one<ConstantType>());
storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>());
return result;
}
template<typename ParametricSparseModelType, typename ConstantType>
ConstantType SamplingModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue() {
invokeSolver();
return this->eqSysResult[this->eqSysInitIndex];
}
template<>
void SamplingModel<storm::models::sparse::Dtmc<storm::RationalFunction>, double>::invokeSolver(){
std::unique_ptr<storm::solver::LinearEquationSolver<double>> solver = storm::utility::solver::LinearEquationSolverFactory<double>().create(this->matrixData.matrix);

7
src/modelchecker/region/SamplingModel.h

@ -61,7 +61,7 @@ namespace storm {
void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
void invokeSolver();
//Some designated states in the model
//Some designated states in the original 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
@ -79,19 +79,18 @@ namespace storm {
* 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
* Finally, there is an assignment that connects every non-constant matrix (or: vector) 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.
*/
std::unordered_map<ParametricType, ConstantType> functions; // the occurring functions together with the corresponding placeholders for the result
struct MatrixData {
storm::storage::SparseMatrix<ConstantType> matrix; //The matrix itself.
std::unordered_map<ParametricType, ConstantType> functions; // the occurring functions together with the corresponding placeholders for the result
std::vector<std::pair<typename storm::storage::SparseMatrix<ConstantType>::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders
} matrixData;
struct VectorData {
std::vector<ConstantType> vector; //The vector itself.
std::unordered_map<ParametricType, ConstantType> functions; // the occurring functions together with the corresponding placeholders for the result
std::vector<std::pair<typename std::vector<ConstantType>::iterator, ConstantType*>> assignment; // Connection of vector entries with placeholders
} vectorData;

Loading…
Cancel
Save