From 67ff27954e33963e32ffe7f6d8c3a2c4e185089c Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 13 Oct 2015 17:27:46 +0200 Subject: [PATCH] refactored approximation model (almost done) Former-commit-id: c7f285906b4daf26de26c50a7dc2ba1b142db880 --- .../region/ApproximationModel.cpp | 525 ++++++++---------- src/modelchecker/region/ApproximationModel.h | 112 ++-- src/modelchecker/region/SamplingModel.cpp | 106 ++-- src/modelchecker/region/SamplingModel.h | 7 +- 4 files changed, 363 insertions(+), 387 deletions(-) diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index d32c99d6f..2f72d61e7 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/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,230 +25,218 @@ namespace storm { namespace region { template - ApproximationModel::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr formula) : formula(formula){ - if(this->formula->isProbabilityOperatorFormula()){ + ApproximationModel::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr 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 probabilityMatrix; - std::vector::index_type> approxRowGroupIndices;// Indicates which rows of the probabilityMatrix belong to a given row of the parametricModel - 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. - //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> rewardModels; - std::vector rewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping - RewTableEntry constantRewEntry; //this value is stored in the rewardEntryToEvalTableMapping for every constant reward vector entry. + STORM_LOG_THROW(parametricModel.hasLabel("target"), storm::exceptions::InvalidArgumentException, "The given Model has no \"target\"-statelabel."); + this->targetStates = parametricModel.getStateLabeling().getStates("target"); + STORM_LOG_THROW(parametricModel.hasLabel("sink"), storm::exceptions::InvalidArgumentException, "The given Model has no \"sink\"-statelabel."); + storm::storage::BitVector sinkStates=parametricModel.getStateLabeling().getStates("sink"); + this->maybeStates = ~(this->targetStates | sinkStates); + STORM_LOG_THROW(parametricModel.getInitialStates().getNumberOfSetBits()==1, storm::exceptions::InvalidArgumentException, "The given model has more or less then one initial state"); + storm::storage::sparse::state_type initialState = *parametricModel.getInitialStates().begin(); + STORM_LOG_THROW(maybeStates.get(initialState), storm::exceptions::InvalidArgumentException, "The value in the initial state of the given model is independent of parameters"); + //The (state-)indices in the equation system will be different from the original ones, as the eq-sys only considers maybestates. + //Hence, we use this vector to translate from old indices to new ones. + std::vector newIndices(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates()); //initialize with some illegal index + std::size_t newIndex=0; + for(auto const& index : maybeStates){ + newIndices[index]=newIndex; + ++newIndex; + } + + //Now pre-compute the information for the equation system. + std::vector rowSubstitutions;// the substitution used in every row + initializeProbabilities(parametricModel, newIndices, rowSubstitutions); if(this->computeRewards){ - std::vector stateActionRewards; - initializeRewards(parametricModel, probabilityMatrix, approxRowGroupIndices, rowSubstitutions, stateActionRewards, rewardEntryToEvalTableMapping, &constantRewEntry); - rewardModels.insert(std::pair>("", storm::models::sparse::StandardRewardModel(boost::optional>(), boost::optional>(std::move(stateActionRewards))))); + initializeRewards(parametricModel, newIndices, rowSubstitutions); } - //Obtain other model ingredients and the approximation model itself - storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); - boost::optional>> noChoiceLabeling; - switch(parametricModel.getType()){ - case storm::models::ModelType::Dtmc: - this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling)); - break; - case storm::models::ModelType::Mdp: - //TODO: use a two-player-game here - this->model=std::make_shared>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling)); - break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Tried to build an Approximation model for an unsupported model type"); + 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(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one() : ConstantType(0.5)); + this->eqSysInitIndex = newIndices[initialState]; + } - //translate the matrixEntryToEvalTableMapping into the actual probability mapping - typename storm::storage::SparseMatrix::index_type matrixRow=0; - auto tableEntry=matrixEntryToEvalTableMapping.begin(); - for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ - for (; matrixRowmodel->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(storm::utility::region::getConstantPart(parEntry.getValue()))); - } else { - this->probabilityMapping.emplace_back(std::make_pair(&((*tableEntry)->second), approxModelEntry)); + template + void ApproximationModel::initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector const& newIndices, std::vector& rowSubstitutions) { + STORM_LOG_DEBUG("Approximation model initialization for probabilities"); + /* 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(); + auto numOfMaybeStates = this->maybeStates.getNumberOfSetBits(); + storm::storage::SparseMatrixBuilder matrixBuilder(numOfMaybeStates, //exact number of rows is unknown at this point, but at least this many + numOfMaybeStates, //columns + 0, //Unknown number of entries + false, // no force dimensions + true, //will have custom row grouping + 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 occurringVariables; + for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ + if(this->maybeStates.get(oldEntry.getColumn())){ + storm::utility::region::gatherOccurringVariables(oldEntry.getValue(), occurringVariables); } - ++approxModelEntry; - ++tableEntry; } - } - } - 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>(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); + std::size_t numOfSubstitutions=1ull< currSubstitution; + std::size_t parameterIndex=0ull; + for(auto const& parameter : occurringVariables){ + if((substitutionId>>parameterIndex)%2==0){ + currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER)); + } + else{ + currSubstitution.insert(std::make_pair(parameter, TypeOfBound::UPPER)); + } + ++parameterIndex; + } + 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 a dummy entry + //Note that this is still executed once, even if no parameters occur. + for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ + if(this->maybeStates.get(oldEntry.getColumn())){ + matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyValue); + } } } + ++curRow; } } - } - - template - void ApproximationModel::initializeProbabilities(ParametricSparseModelType const& parametricModel, - storm::storage::SparseMatrix& probabilityMatrix, - std::vector::index_type>& approxRowGroupIndices, - std::vector& rowSubstitutions, - std::vector& matrixEntryToEvalTableMapping, - ProbTableEntry* constantEntry) { - 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 matrixBuilder(0, //unknown number of rows - parametricModel.getTransitionMatrix().getColumnCount(), - 0, //Unknown number of entries - true, //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::index_type apprModelMatrixRow=0; - for(typename storm::storage::SparseMatrix::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); - } - // 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(paramRow)){ - storm::utility::region::gatherOccurringVariables(entry.getValue(), occurringVariables); - } - std::size_t numOfSubstitutions=1ull< currSubstitution; - std::size_t parameterIndex=0; - for(auto const& parameter : occurringVariables){ - if((substitutionId>>parameterIndex)%2==0){ - currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER)); - } - else{ - currSubstitution.insert(std::make_pair(parameter, TypeOfBound::UPPER)); + 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(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(storm::utility::zero()); + if(!this->computeRewards){ + for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ + if(this->targetStates.get(oldEntry.getColumn())){ + targetProbability += oldEntry.getValue(); + } } - ++parameterIndex; } - std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->probabilitySubstitutions, 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 - //Note that this is still executed once, even if no parameters occur. - ConstantType dummyEntry=storm::utility::one(); - for(auto const& entry : parametricModel.getTransitionMatrix().getRow(paramRow)){ - if(storm::utility::isConstant(entry.getValue())){ - matrixEntryToEvalTableMapping.emplace_back(constantEntry); - } else { - ++numOfNonConstEntries; - auto evalTableIt = this->probabilityEvaluationTable.insert(ProbTableEntry(FunctionSubstitution(entry.getValue(), substitutionIndex), storm::utility::zero())).first; - matrixEntryToEvalTableMapping.emplace_back(&(*evalTableIt)); + //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(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; + } } - matrixBuilder.addNextValue(apprModelMatrixRow, entry.getColumn(), dummyEntry); - dummyEntry=storm::utility::zero(); + if(!this->computeRewards){ + if(storm::utility::isConstant(targetProbability)){ + *vectorIt = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability)); + } else { + 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; + } + } + ++vectorIt; } - ++apprModelMatrixRow; + ++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 - void ApproximationModel::initializeRewards(ParametricSparseModelType const& parametricModel, - storm::storage::SparseMatrix const& probabilityMatrix, - std::vector::index_type> const& approxRowGroupIndices, - std::vector const& rowSubstitutions, - std::vector& stateActionRewardVector, - std::vector& rewardEntryToEvalTableMapping, - RewTableEntry* constantEntry){ + void ApproximationModel::initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices, std::vector 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 occurringRewVariables; - std::set occurringProbVariables; - if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])){ - ConstantType reward = storm::utility::region::convertNumber(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(); + 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(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]; matrixRowmatrixData.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 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::iterator functionsIt; + for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRowmatrixData.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 rewardSubstitution; + //Get the correct substitution for this matrixRow + std::map substitution; for(auto const& rewardVar : occurringRewVariables){ - 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)); + typename std::map::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(storm::utility::zero(), storm::utility::zero())) - ).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()); - 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,16 +247,32 @@ namespace storm { } template - std::shared_ptr> const& ApproximationModel::getModel() const { - return this->model; + std::vector ApproximationModel::computeValues(ParameterRegion const& region, bool computeLowerBounds) { + instantiate(region, computeLowerBounds); + invokeSolver(computeLowerBounds); + + std::vector result(this->maybeStates.size()); + storm::utility::vector::setVectorValues(result, this->maybeStates, this->eqSysResult); + storm::utility::vector::setVectorValues(result, this->targetStates, this->computeRewards ? storm::utility::zero() : storm::utility::one()); + storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity() : storm::utility::zero()); + + return result; } template - void ApproximationModel::instantiate(const ParameterRegion& region) { + ConstantType ApproximationModel::computeInitialStateValue(ParameterRegion const& region, bool computeLowerBounds) { + instantiate(region, computeLowerBounds); + invokeSolver(computeLowerBounds); + return this->eqSysResult[this->eqSysInitIndex]; + } + + template + void ApproximationModel::instantiate(const ParameterRegion& region, bool computeLowerBounds) { //Instantiate the substitutions - std::vector> instantiatedSubs(this->probabilitySubstitutions.size()); - for(uint_fast64_t substitutionIndex=0; substitutionIndexprobabilitySubstitutions.size(); ++substitutionIndex){ - for(std::pair const& sub : this->probabilitySubstitutions[substitutionIndex]){ + std::vector> instantiatedSubs(this->funcSubData.substitutions.size()); + std::vector> choseOptimalParameters(this->funcSubData.substitutions.size()); + for(std::size_t substitutionIndex=0; substitutionIndexfuncSubData.substitutions.size(); ++substitutionIndex){ + for(std::pair const& sub : this->funcSubData.substitutions[substitutionIndex]){ switch(sub.second){ case TypeOfBound::LOWER: instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first))); @@ -275,117 +280,65 @@ namespace storm { case TypeOfBound::UPPER: instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getUpperBound(sub.first))); break; + case TypeOfBound::CHOSEOPTIMAL: + //Insert some dummy value + instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, storm::utility::one())); + choseOptimalParameters[substitutionIndex].insert(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( - 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)); - } - - 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 function+substitution results into placeholders + for(auto& functionResult : this->funcSubData.functions){ + auto& funcSub = functionResult.first; + auto& result = functionResult.second; + result = computeLowerBounds ? storm::utility::infinity() : storm::utility::zero(); + //Iterate over the different combinations of lower and upper bounds and update the min and max values + auto const& vertices=region.getVerticesOfRegion(this->choseOptimalParameters[funcSub.getSubstitution()]); + for(auto const& vertex : vertices){ + //extend the substitution + for(auto const& vertexSub : vertex){ + instantiatedSubs[funcSub.getSubstitution()][vertexSub.first]=vertexSub.second; } - } - //write entries into evaluation table - for(auto& tableEntry : this->rewardEvaluationTable){ - auto& funcSub = tableEntry.first; - auto& minMax = tableEntry.second; - minMax.first=storm::utility::infinity(); - minMax.second=storm::utility::zero(); - //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()]); - for(auto const& vertex : vertices){ - //extend the substitution - for(auto const& vertexSub : vertex){ - instantiatedRewardSubs[funcSub.getSubstitution()][vertexSub.first]=vertexSub.second; - } - ConstantType currValue = storm::utility::region::convertNumber( - storm::utility::region::evaluateFunction( - funcSub.getFunction(), - instantiatedRewardSubs[funcSub.getSubstitution()] + //evaluate the function + ConstantType currValue = storm::utility::region::convertNumber( + storm::utility::region::evaluateFunction( + funcSub.getFunction(), + instantiatedSubs[funcSub.getSubstitution()] ) ); - minMax.first=std::min(minMax.first, currValue); - minMax.second=std::max(minMax.second, currValue); - } + result = computeLowerBounds ? std::min(result, currValue) : std::max(result, currValue); } - //Note: the rewards are written to the model as soon as the optimality type is known (i.e. in computeValues) - } - } - - template - std::unique_ptr ApproximationModel::computeValues(storm::solver::OptimizationDirection const& approximationOpDir) { - std::unique_ptr modelChecker; - switch(this->getModel()->getType()){ - case storm::models::ModelType::Mdp: - modelChecker = std::make_unique>>(*this->model->template as>()); - 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"); } - std::unique_ptr resultPtr; - //TODO: use this when two player games are available - boost::optional 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); - } - 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 noRewardModelName; //it should be uniquely given - resultPtr = modelChecker->computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula(), noRewardModelName, false, approximationOpDir); + + //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)); } - else { - //perform model checking - resultPtr = modelChecker->computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula(), false, approximationOpDir); + for(auto& assignment : this->vectorData.assignment){ + *assignment.first=*assignment.second; } - return resultPtr; } + + + template<> + void ApproximationModel, double>::invokeSolver(bool computeLowerBounds){ + storm::solver::SolveGoal goal(computeLowerBounds); + std::unique_ptr> solver = storm::solver::configureMinMaxLinearEquationSolver(goal, storm::utility::solver::MinMaxLinearEquationSolverFactory(), this->matrixData.matrix); + solver->solveEquationSystem(this->eqSysResult, this->vectorData.vector); + } + + template<> + void ApproximationModel, double>::invokeSolver(bool computeLowerBounds){ + storm::solver::SolveGoal player2Goal(computeLowerBounds); + std::unique_ptr> solver = storm::utility::solver::GameSolverFactory().create(this->player1Matrix, this->matrixData.matrix); + solver->solveGame(this->player1Goal.direction(), player2Goal.direction(), this->eqSysResult, this->vectorData.vector); + } + + #ifdef STORM_HAVE_CARL diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h index c79c50bc3..4b11f22f2 100644 --- a/src/modelchecker/region/ApproximationModel.h +++ b/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 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> const& getModel() const; + std::vector computeValues(ParameterRegion 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 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 computeValues(storm::solver::OptimizationDirection const& approximationOpDir); - + ConstantType computeInitialStateValue(ParameterRegion const& region, bool computeLowerBounds); private: //This enum helps to store how a parameter will be substituted. @@ -121,43 +122,60 @@ namespace storm { } }; - typedef typename std::unordered_map::value_type ProbTableEntry; - typedef typename std::unordered_map, FuncSubHash>::value_type RewTableEntry; - - void initializeProbabilities(ParametricSparseModelType const& parametricModel, storm::storage::SparseMatrix& probabilityMatrix, std::vector::index_type>& approxRowGroupIndices, std::vector& rowSubstitutions, std::vector& matrixEntryToEvalTableMapping, ProbTableEntry* constantEntry); - void initializeRewards(ParametricSparseModelType const& parametricModel, storm::storage::SparseMatrix const& probabilityMatrix, std::vector::index_type> const& approxRowGroupIndices, std::vector const& rowSubstitutions, std::vector& stateActionRewardVector, std::vector& rewardEntryToEvalTableMapping, RewTableEntry* constantEntry); - - //The Model with which we work - std::shared_ptr> model; - //The formula for which we will compute the values - std::shared_ptr formula; + typedef typename std::unordered_map::value_type FunctionEntry; + // typedef typename std::unordered_map, FuncSubHash>::value_type RewTableEntry; + + void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector const& newIndices, std::vector& rowSubstitutions); + void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices, std::vector const& rowSubstitutions); + void initializePlayer1Matrix(ParametricSparseModelType const& parametricModel, std::vector const& newIndices); + void instantiate(ParameterRegion 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 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; - - // 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 probabilityEvaluationTable; - //For rewards, we map to the minimal value and the maximal value (depending on the CHOSEOPTIMAL parameters). - std::unordered_map, FuncSubHash> rewardEvaluationTable; - - //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region) - 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 (i.e. CHOSEOPTIMAL parameters) - std::vector> rewardSubstitutions; - std::vector> 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::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::iterator>> rewardMapping; - + //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 const& player1Matrix; + + /* 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 functions; + // the occurring reward functions together with the corresponding placeholders for the minimal and maximal result + // //std::unordered_map, FuncSubHash> rewFunctions; + //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region) + std::vector> 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> rewSubs; + // //std::vector> choseOptimalRewardPars; + } funcSubData; + struct MatrixData { + storm::storage::SparseMatrix matrix; //The matrix itself. + std::vector::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders + // std::vector::iterator>> rewardMapping; + } matrixData; + struct VectorData { + std::vector vector; //The vector itself. + std::vector::iterator, ConstantType*>> assignment; // Connection of vector entries with placeholders + } vectorData; + }; } //namespace region } diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index 8e3b03107..1b091b440 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/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(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one() : ConstantType(0.5)); this->eqSysInitIndex = newIndices[initialState]; - this->solveGoal = storm::solver::SolveGoal(storm::logic::isLowerBound(formula->getComparisonType())); } template void SamplingModel::initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector 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(); 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 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(storm::utility::zero()); for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ if(this->maybeStates.get(oldEntry.getColumn())){ //check if we need to insert a diagonal entry @@ -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()); } - if(!this->computeRewards){ - if(storm::utility::isConstant(targetProbability)){ - this->vectorData.vector.push_back(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability))); - } else { - typename std::unordered_map::iterator functionsIt = this->vectorData.functions.insert(FunctionEntry(targetProbability, dummyValue)).first; - this->vectorData.assignment.emplace_back(std::make_pair(this->vectorData.vector.end(), &(functionsIt->second))); - this->vectorData.vector.push_back(dummyValue); - } - } ++curRow; } } 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(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(storm::utility::zero()); 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::iterator functionsIt; if(isDtmc){ if(eqSysMatrixEntry->getColumn()==newIndices[oldRowGroup]){ //Diagonal entries get 1-f(x) - functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::one()-oldEntry.getValue(), dummyValue)).first; + functionsIt = this->functions.insert(FunctionEntry(storm::utility::one()-oldEntry.getValue(), dummyValue)).first; } else { - functionsIt = this->matrixData.functions.insert(FunctionEntry(storm::utility::zero()-oldEntry.getValue(), dummyValue)).first; + functionsIt = this->functions.insert(FunctionEntry(storm::utility::zero()-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(storm::utility::region::getConstantPart(targetProbability)); + } else { + typename std::unordered_map::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 void SamplingModel::initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices){ - //Run through the state reward vector... - this->vectorData.vector.reserve(this->maybeStates.getNumberOfSetBits()); //Important! iterators have to remain valid - this->vectorData.assignment.reserve(this->maybeStates.getNumberOfSetBits()); + //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(); + 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(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state]))); + *vectorIt = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state])); } else { - typename std::unordered_map::iterator functionsIt = this->vectorData.functions.insert(FunctionEntry(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], dummyValue)).first; - this->vectorData.assignment.emplace_back(std::make_pair(this->vectorData.vector.end(), &(functionsIt->second))); - this->vectorData.vector.push_back(dummyValue); + typename std::unordered_map::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 @@ -198,14 +204,31 @@ namespace storm { //Intentionally left empty } + template + std::vector SamplingModel::computeValues() { + invokeSolver(); + std::vector result(this->maybeStates.size()); + storm::utility::vector::setVectorValues(result, this->maybeStates, this->eqSysResult); + storm::utility::vector::setVectorValues(result, this->targetStates, this->computeRewards ? storm::utility::zero() : storm::utility::one()); + storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity() : storm::utility::zero()); + + return result; + } + + template + ConstantType SamplingModel::computeInitialStateValue() { + invokeSolver(); + return this->eqSysResult[this->eqSysInitIndex]; + } + template void SamplingModel::instantiate(std::mapconst& point) { //write results into the placeholders - for(auto& functionResult : this->matrixData.functions){ + for(auto& functionResult : this->functions){ functionResult.second=storm::utility::region::convertNumber( storm::utility::region::evaluateFunction(functionResult.first, point)); } - for(auto& functionResult : this->vectorData.functions){ + for(auto& functionResult : this->functions){ functionResult.second=storm::utility::region::convertNumber( storm::utility::region::evaluateFunction(functionResult.first, point)); } @@ -218,23 +241,6 @@ namespace storm { *assignment.first=*assignment.second; } } - - template - std::vector SamplingModel::computeValues() { - invokeSolver(); - std::vector result(this->maybeStates.size()); - storm::utility::vector::setVectorValues(result, this->maybeStates, this->eqSysResult); - storm::utility::vector::setVectorValues(result, this->targetStates, this->computeRewards ? storm::utility::zero() : storm::utility::one()); - storm::utility::vector::setVectorValues(result, ~(this->maybeStates | this->targetStates), this->computeRewards ? storm::utility::infinity() : storm::utility::zero()); - - return result; - } - - template - ConstantType SamplingModel::computeInitialStateValue() { - invokeSolver(); - return this->eqSysResult[this->eqSysInitIndex]; - } template<> void SamplingModel, double>::invokeSolver(){ diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index 3fec3b086..42b8bb361 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -61,7 +61,7 @@ namespace storm { void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector 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 functions; // the occurring functions together with the corresponding placeholders for the result struct MatrixData { storm::storage::SparseMatrix matrix; //The matrix itself. - std::unordered_map functions; // the occurring functions together with the corresponding placeholders for the result std::vector::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders } matrixData; struct VectorData { std::vector vector; //The vector itself. - std::unordered_map functions; // the occurring functions together with the corresponding placeholders for the result std::vector::iterator, ConstantType*>> assignment; // Connection of vector entries with placeholders } vectorData;