Browse Source

more efficient instantiation of matrices.

Former-commit-id: 0f4c36f2f2
tempestpy_adaptions
TimQu 9 years ago
parent
commit
28326b14e4
  1. 458
      src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp
  2. 69
      src/modelchecker/reachability/SparseDtmcRegionModelChecker.h
  3. 170
      src/modelchecker/region/ApproximationModel.cpp
  4. 82
      src/modelchecker/region/ApproximationModel.h
  5. 107
      src/modelchecker/region/SamplingModel.cpp
  6. 72
      src/modelchecker/region/SamplingModel.h
  7. 5
      src/utility/regions.cpp
  8. 10
      src/utility/regions.h

458
src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp

@ -4,7 +4,6 @@
#include "src/adapters/CarlAdapter.h"
//#include "src/storage/StronglyConnectedComponentDecomposition.h"
#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@ -15,7 +14,9 @@
#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h"
//#include "modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
#include "src/modelchecker/region/ApproximationModel.h"
#include "src/modelchecker/region/SamplingModel.h"
#include "src/exceptions/InvalidPropertyException.h"
#include "src/exceptions/InvalidStateException.h"
@ -183,9 +184,10 @@ namespace storm {
eliminationModelChecker(model),
smtSolver(nullptr),
probabilityOperatorFormula(nullptr),
sampleDtmc(nullptr),
approxMdp(nullptr),
isReachProbFunctionComputed(false) {
samplingModel(nullptr),
approximationModel(nullptr),
isReachProbFunctionComputed(false),
isResultConstant(false) {
//intentionally left empty
}
@ -221,66 +223,19 @@ namespace storm {
std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now();
STORM_LOG_THROW(this->canHandle(formula), storm::exceptions::IllegalArgumentException, "Tried to specify a formula that can not be handled.");
//Get subformula, initial state, target states
//Get subformula, target states
//Note: canHandle already ensures that the formula has the right shape and that the model has a single initial state.
this->probabilityOperatorFormula= std::unique_ptr<storm::logic::ProbabilityOperatorFormula>(new storm::logic::ProbabilityOperatorFormula(formula.asStateFormula().asProbabilityOperatorFormula()));
storm::logic::EventuallyFormula const& eventuallyFormula = this->probabilityOperatorFormula->getSubformula().asPathFormula().asEventuallyFormula();
std::unique_ptr<CheckResult> targetStatesResultPtr = this->eliminationModelChecker.check(eventuallyFormula.getSubformula());
storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector();
// Then, compute the subset of states that has a probability of 0 or 1, respectively.
std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(model, storm::storage::BitVector(model.getNumberOfStates(),true), targetStates);
storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
// If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction.
if (model.getInitialStates().isDisjointFrom(maybeStates)) {
STORM_LOG_WARN("The probability of the initial state is constant (0 or 1)");
this->reachProbFunction = statesWithProbability0.get(*model.getInitialStates().begin()) ? storm::utility::zero<ParametricType>() : storm::utility::one<ParametricType>();
this->isReachProbFunctionComputed=true;
}
// Determine the set of states that is reachable from the initial state without jumping over a target state.
storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, statesWithProbability1);
// Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
maybeStates &= reachableStates;
// Create a vector for the probabilities to go to a state with probability 1 in one step.
this->oneStepProbabilities = model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
// Determine the initial state of the sub-model.
//storm::storage::BitVector newInitialStates = model.getInitialStates() % maybeStates;
this->initialState = *(model.getInitialStates() % maybeStates).begin();
// We then build the submatrix that only has the transitions of the maybe states.
storm::storage::SparseMatrix<ParametricType> submatrix = model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
storm::storage::SparseMatrix<ParametricType> submatrixTransposed = submatrix.transpose();
// Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
this->flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix);
this->flexibleBackwardTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrixTransposed, true);
// Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated.
this->subsystem = storm::storage::BitVector(submatrix.getRowCount(), true);
std::chrono::high_resolution_clock::time_point timeInitialStateEliminationStart = std::chrono::high_resolution_clock::now();
// eliminate all states with only constant outgoing transitions
//TODO: maybe also states with constant incoming tranistions. THEN the ordering of the eliminated states does matter.
eliminateStatesConstSucc(this->subsystem, this->flexibleTransitions, this->flexibleBackwardTransitions, this->oneStepProbabilities, this->hasOnlyLinearFunctions, this->initialState);
STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl);
std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl;
//eliminate the remaining states to get the reachability probability function
this->sparseTransitions = this->flexibleTransitions.getSparseMatrix();
this->sparseBackwardTransitions = this->sparseTransitions.transpose();
std::chrono::high_resolution_clock::time_point timeInitialStateEliminationEnd = std::chrono::high_resolution_clock::now();
initializeSampleDtmcAndApproxMdp(this->sampleDtmc,this->sampleDtmcMapping,this->approxMdp,this->approxMdpMapping,this->approxMdpSubstitutions,this->subsystem, this->sparseTransitions,this->oneStepProbabilities, this->initialState);
computeSimplifiedModel(targetStates);
initializeSampleAndApproxModel();
//some information for statistics...
std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now();
this->timePreprocessing= timePreprocessingEnd - timePreprocessingStart;
this->timeInitialStateElimination = timeInitialStateEliminationEnd-timeInitialStateEliminationStart;
this->numOfCheckedRegions=0;
this->numOfRegionsSolvedThroughSampling=0;
this->numOfRegionsSolvedThroughApproximation=0;
@ -298,42 +253,134 @@ namespace storm {
}
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::eliminateStatesConstSucc(
storm::storage::BitVector& subsys,
FlexibleMatrix& flexTransitions,
FlexibleMatrix& flexBackwardTransitions,
std::vector<ParametricType>& oneStepProbs,
bool& allFunctionsAreLinear,
storm::storage::sparse::state_type const& initialState) {
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::computeSimplifiedModel(storm::storage::BitVector const& targetStates){
//Compute the subset of states that have a probability of 0 or 1, respectively and reduce the considered states accordingly.
std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->model, storm::storage::BitVector(this->model.getNumberOfStates(),true), targetStates);
storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
// If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction.
if (this->model.getInitialStates().isDisjointFrom(maybeStates)) {
STORM_LOG_WARN("The probability of the initial state is constant (0 or 1)");
this->reachProbFunction = statesWithProbability0.get(*(this->model.getInitialStates().begin())) ? storm::utility::zero<ParametricType>() : storm::utility::one<ParametricType>();
this->isReachProbFunctionComputed=true;
this->isResultConstant=true;
}
// Determine the set of states that is reachable from the initial state without jumping over a target state.
storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->model.getTransitionMatrix(), this->model.getInitialStates(), maybeStates, statesWithProbability1);
// Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
maybeStates &= reachableStates;
// Create a vector for the probabilities to go to a state with probability 1 in one step.
std::vector<ParametricType> oneStepProbabilities = this->model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
// Determine the initial state of the sub-model.
storm::storage::sparse::state_type initialState = *(this->model.getInitialStates() % maybeStates).begin();
// We then build the submatrix that only has the transitions of the maybe states.
storm::storage::SparseMatrix<ParametricType> submatrix = this->model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
// eliminate all states with only constant outgoing transitions
//TODO: maybe also states with constant incoming tranistions. THEN the ordering of the eliminated states does matter.
std::chrono::high_resolution_clock::time_point timeInitialStateEliminationStart = std::chrono::high_resolution_clock::now();
// Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
FlexibleMatrix flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix);
FlexibleMatrix flexibleBackwardTransitions= this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix.transpose(), true);
// Create a bit vector that represents the current subsystem, i.e., states that we have not eliminated.
storm::storage::BitVector subsystem = storm::storage::BitVector(submatrix.getRowCount(), true);
//temporarily unselect the initial state to skip it.
subsys.set(initialState, false);
allFunctionsAreLinear=true;
subsystem.set(initialState, false);
this->hasOnlyLinearFunctions=true;
boost::optional<std::vector<ParametricType>> missingStateRewards;
for (auto const& state : subsys) {
for (auto const& state : subsystem) {
bool onlyConstantOutgoingTransitions=true;
for(auto& entry : flexTransitions.getRow(state)){
for(auto const& entry : submatrix.getRow(state)){
if(!this->parametricTypeComparator.isConstant(entry.getValue())){
onlyConstantOutgoingTransitions=false;
allFunctionsAreLinear &= storm::utility::regions::functionIsLinear<ParametricType>(entry.getValue());
this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(entry.getValue());
}
}
if(!this->parametricTypeComparator.isConstant(oneStepProbs[state])){
if(!this->parametricTypeComparator.isConstant(oneStepProbabilities[state])){
onlyConstantOutgoingTransitions=false;
allFunctionsAreLinear &= storm::utility::regions::functionIsLinear<ParametricType>(oneStepProbs[state]);
this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear<ParametricType>(oneStepProbabilities[state]);
}
if(onlyConstantOutgoingTransitions){
this->eliminationModelChecker.eliminateState(flexTransitions, oneStepProbs, state, flexBackwardTransitions, missingStateRewards);
subsys.set(state,false);
this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards);
subsystem.set(state,false);
}
}
subsys.set(initialState, true);
subsystem.set(initialState, true);
std::chrono::high_resolution_clock::time_point timeInitialStateEliminationEnd = std::chrono::high_resolution_clock::now();
this->timeInitialStateElimination = timeInitialStateEliminationEnd-timeInitialStateEliminationStart;
STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl);
std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl;
//Build the simplified model
//The matrix. The flexibleTransitions matrix might have empty rows where states have been eliminated.
//The new matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
std::vector<storm::storage::sparse::state_type> newStateIndexMap(flexibleTransitions.getNumberOfRows(), flexibleTransitions.getNumberOfRows()); //initialize with some illegal index to easily check if a transition leads to an unselected state
storm::storage::sparse::state_type newStateIndex=0;
for(auto const& state : subsystem){
newStateIndexMap[state]=newStateIndex;
++newStateIndex;
}
//We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead
storm::storage::sparse::state_type numStates=newStateIndex+2;
storm::storage::sparse::state_type targetState=numStates-2;
storm::storage::sparse::state_type sinkState= numStates-1;
//We can now fill in the data.
storm::storage::SparseMatrixBuilder<ParametricType> matrixBuilder(numStates, numStates);
for(storm::storage::sparse::state_type oldStateIndex : subsystem){
ParametricType valueToSinkState=storm::utility::regions::getNewFunction<ParametricType, CoefficientType>(storm::utility::one<CoefficientType>());
//go through columns:
for(auto const& entry: flexibleTransitions.getRow(oldStateIndex)){
STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=flexibleTransitions.getNumberOfRows(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated.");
valueToSinkState-=entry.getValue();
matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],entry.getValue());
}
//transition to target state
if(!this->parametricTypeComparator.isZero(oneStepProbabilities[oldStateIndex])){
valueToSinkState-=oneStepProbabilities[oldStateIndex];
matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, oneStepProbabilities[oldStateIndex]);
}
//transition to sink state
if(!this->parametricTypeComparator.isZero(valueToSinkState)){
matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, valueToSinkState);
}
}
//add self loops on the additional states (i.e., target and sink)
matrixBuilder.addNextValue(targetState, targetState, storm::utility::one<ParametricType>());
matrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one<ParametricType>());
//The labeling
storm::models::sparse::StateLabeling labeling(numStates);
storm::storage::BitVector initLabel(numStates, false);
initLabel.set(newStateIndexMap[initialState], true);
labeling.addLabel("init", std::move(initLabel));
storm::storage::BitVector targetLabel(numStates, false);
targetLabel.set(targetState, true);
labeling.addLabel("target", std::move(targetLabel));
storm::storage::BitVector sinkLabel(numStates, false);
sinkLabel.set(sinkState, true);
labeling.addLabel("sink", std::move(sinkLabel));
// other ingredients
boost::optional<std::vector<ParametricType>> noStateRewards;
boost::optional<storm::storage::SparseMatrix<ParametricType>> noTransitionRewards;
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;
// the final model
this->simplifiedModel = std::make_shared<storm::models::sparse::Dtmc<ParametricType>>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
}
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSampleDtmcAndApproxMdp(
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSampleAndApproxModel(){
//now create the models
this->approximationModel=std::make_shared<ApproximationModel>(*this->simplifiedModel);
this->samplingModel=std::make_shared<SamplingModel>(*this->simplifiedModel);
}
/*
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSampleDtmcAndApproxMdp2(
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>>& sampleDtmc,
std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>> & sampleDtmcMapping,
std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> & approxMdp,
@ -344,7 +391,6 @@ namespace storm {
std::vector<ParametricType> const& oneStepProbs,
storm::storage::sparse::state_type const& initState) {
std::set<ParametricType> functionSet;
//the matrix "transitions" might have empty rows where states have been eliminated.
//The DTMC/ MDP matrices should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
@ -378,7 +424,6 @@ namespace storm {
for(auto const& entry: transitions.getRow(oldStateIndex)){
valueToSinkState-=entry.getValue();
storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringParameters);
functionSet.insert(entry.getValue());
dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],dummyEntry);
mdpMatrixBuilder.addNextValue(currentMdpRow, newStateIndexMap[entry.getColumn()], dummyEntry);
rowInformation.emplace_back(std::make_pair(newStateIndexMap[entry.getColumn()], dummyEntry));
@ -389,7 +434,6 @@ namespace storm {
if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){
valueToSinkState-=oneStepProbs[oldStateIndex];
storm::utility::regions::gatherOccurringVariables(oneStepProbs[oldStateIndex], occurringParameters);
functionSet.insert(oneStepProbs[oldStateIndex]);
dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, dummyEntry);
mdpMatrixBuilder.addNextValue(currentMdpRow, targetState, dummyEntry);
rowInformation.emplace_back(std::make_pair(targetState, dummyEntry));
@ -401,7 +445,6 @@ namespace storm {
mdpMatrixBuilder.addNextValue(currentMdpRow, sinkState, dummyEntry);
rowInformation.emplace_back(std::make_pair(sinkState, dummyEntry));
oneStepToSink[oldStateIndex]=valueToSinkState;
functionSet.insert(valueToSinkState);
}
// Find the different combinations of occuring variables and lower/upper bounds (mathematically speaking we want to have the set 'occuringVariables times {u,l}' )
std::size_t numOfSubstitutions=std::pow(2,occurringParameters.size()); //1 substitution = 1 combination of lower and upper bounds
@ -457,7 +500,6 @@ namespace storm {
// the final dtmc and mdp
sampleDtmc=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(dtmcMatrixBuilder.build(), std::move(dtmcStateLabeling), std::move(dtmcNoStateRewards), std::move(dtmcNoTransitionRewards), std::move(dtmcNoChoiceLabeling));
approxMdp=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(mdpMatrixBuilder.build(), std::move(mdpStateLabeling), std::move(mdpNoStateRewards), std::move(mdpNoTransitionRewards), std::move(mdpNoChoiceLabeling));
std::cout << "THERE ARE " << functionSet.size() << " DIFFERENT RATIONAL FUNCTIONS" << std::endl;
//build mapping vectors and the substitution vector.
sampleDtmcMapping=std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>>();
sampleDtmcMapping.reserve(sampleDtmc->getTransitionMatrix().getEntryCount());
@ -532,12 +574,29 @@ namespace storm {
}
}
}
*/
template<typename ParametricType, typename ConstantType>
ParametricType SparseDtmcRegionModelChecker<ParametricType, ConstantType>::getReachProbFunction() {
if(!this->isReachProbFunctionComputed){
std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionStart = std::chrono::high_resolution_clock::now();
this->reachProbFunction = computeReachProbFunction(this->subsystem, this->flexibleTransitions, this->flexibleBackwardTransitions, this->sparseTransitions, this->sparseBackwardTransitions, this->oneStepProbabilities, this->initialState);
//get the one step probabilities and the transition matrix of the simplified model without target/sink state
//TODO check if this takes long... we could also store the oneSteps while creating the simplified model. Or(?) we keep the matrix the way it is and give the target state one step probability 1.
storm::storage::SparseMatrix<ParametricType> backwardTransitions(this->simplifiedModel->getBackwardTransitions());
std::vector<ParametricType> oneStepProbabilities(this->simplifiedModel->getNumberOfStates()-2, storm::utility::zero<ParametricType>());
for(auto const& entry : backwardTransitions.getRow(*(this->simplifiedModel->getStates("target").begin()))){
if(entry.getColumn()<oneStepProbabilities.size()){
oneStepProbabilities[entry.getColumn()]=entry.getValue();
} //else case: only holds for the entry that corresponds to the selfloop on the target state..
}
storm::storage::BitVector maybeStates=~(this->simplifiedModel->getStates("target") | this->simplifiedModel->getStates("sink"));
backwardTransitions=backwardTransitions.getSubmatrix(false,maybeStates,maybeStates);
storm::storage::SparseMatrix<ParametricType> forwardTransitions=this->simplifiedModel->getTransitionMatrix().getSubmatrix(false,maybeStates,maybeStates);
//now compute the functions using methods from elimination model checker
storm::storage::BitVector newInitialStates = this->simplifiedModel->getInitialStates() % maybeStates;
storm::storage::BitVector phiStates(this->simplifiedModel->getNumberOfStates(), true);
boost::optional<std::vector<ParametricType>> noStateRewards;
std::vector<std::size_t> statePriorities = this->eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities);
this->reachProbFunction=this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, this->simplifiedModel->getStates("target"), noStateRewards, statePriorities);
std::chrono::high_resolution_clock::time_point timeComputeReachProbFunctionEnd = std::chrono::high_resolution_clock::now();
std::string funcStr = " (/ " +
this->reachProbFunction.nominator().toString(false, true) + " " +
@ -551,64 +610,7 @@ namespace storm {
return this->reachProbFunction;
}
template<typename ParametricType, typename ConstantType>
ParametricType SparseDtmcRegionModelChecker<ParametricType, ConstantType>::computeReachProbFunction(
storm::storage::BitVector const& subsys,
FlexibleMatrix const& flexTransitions,
FlexibleMatrix const& flexBackwardTransitions,
storm::storage::SparseMatrix<ParametricType> const& spTransitions,
storm::storage::SparseMatrix<ParametricType> const& spBackwardTransitions,
std::vector<ParametricType> const& oneStepProbs,
storm::storage::sparse::state_type const& initState){
//Note: this function is very similar to eliminationModelChecker.computeReachabilityValue
//get working copies of the flexible transitions and oneStepProbabilities with which we will actually work (eliminate states etc).
FlexibleMatrix workingCopyFlexTrans(flexTransitions);
FlexibleMatrix workingCopyFlexBackwTrans(flexBackwardTransitions);
std::vector<ParametricType> workingCopyOneStepProbs(oneStepProbs);
storm::storage::BitVector initialStates(subsys.size(),false);
initialStates.set(initState, true);
std::vector<std::size_t> statePriorities = this->eliminationModelChecker.getStatePriorities(spTransitions, spBackwardTransitions, initialStates, workingCopyOneStepProbs);
boost::optional<std::vector<ParametricType>> missingStateRewards;
// Remove the initial state from the states which we need to eliminate.
storm::storage::BitVector statesToEliminate(subsys);
statesToEliminate.set(initState, false);
//pure state elimination or hybrid method?
if (storm::settings::sparseDtmcEliminationModelCheckerSettings().getEliminationMethod() == storm::settings::modules::SparseDtmcEliminationModelCheckerSettings::EliminationMethod::State) {
std::vector<storm::storage::sparse::state_type> states(statesToEliminate.begin(), statesToEliminate.end());
std::sort(states.begin(), states.end(), [&statePriorities] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statePriorities[a] < statePriorities[b]; });
STORM_LOG_DEBUG("Eliminating " << states.size() << " states using the state elimination technique." << std::endl);
for (auto const& state : states) {
this->eliminationModelChecker.eliminateState(workingCopyFlexTrans, workingCopyOneStepProbs, state, workingCopyFlexBackwTrans, missingStateRewards);
}
}
else if (storm::settings::sparseDtmcEliminationModelCheckerSettings().getEliminationMethod() == storm::settings::modules::SparseDtmcEliminationModelCheckerSettings::EliminationMethod::Hybrid) {
uint_fast64_t maximalDepth = 0;
std::vector<storm::storage::sparse::state_type> entryStateQueue;
STORM_LOG_DEBUG("Eliminating " << statesToEliminate.size() << " states using the hybrid elimination technique." << std::endl);
maximalDepth = eliminationModelChecker.treatScc(workingCopyFlexTrans, workingCopyOneStepProbs, initialStates, statesToEliminate, spTransitions, workingCopyFlexBackwTrans, false, 0, storm::settings::sparseDtmcEliminationModelCheckerSettings().getMaximalSccSize(), entryStateQueue, missingStateRewards, statePriorities);
// If the entry states were to be eliminated last, we need to do so now.
STORM_LOG_DEBUG("Eliminating " << entryStateQueue.size() << " entry states as a last step.");
if (storm::settings::sparseDtmcEliminationModelCheckerSettings().isEliminateEntryStatesLastSet()) {
for (auto const& state : entryStateQueue) {
eliminationModelChecker.eliminateState(workingCopyFlexTrans, workingCopyOneStepProbs, state, workingCopyFlexBackwTrans, missingStateRewards);
}
}
}
else {
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "The selected state elimination Method has not been implemented.");
}
//finally, eliminate the initial state
this->eliminationModelChecker.eliminateState(workingCopyFlexTrans, workingCopyOneStepProbs, initState, workingCopyFlexBackwTrans, missingStateRewards);
return workingCopyOneStepProbs[initState];
}
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, const ParametricType& reachProbFunc, const storm::logic::ProbabilityOperatorFormula& formula) {
@ -788,18 +790,8 @@ namespace storm {
valueInBoundOfFormula = this->valueIsInBoundOfFormula(storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(getReachProbFunction(), point));
}
else{
//put the values into the dtmc matrix
for( std::pair<ParametricType, typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&>& mappingPair : this->sampleDtmcMapping){
mappingPair.second.setValue(storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(mappingPair.first, point)
)
);
}
std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr);
storm::modelchecker::SparseDtmcPrctlModelChecker<ConstantType> modelChecker(*this->sampleDtmc);
std::unique_ptr<storm::modelchecker::CheckResult> resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false);
valueInBoundOfFormula=this->valueIsInBoundOfFormula(resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector()[*this->sampleDtmc->getInitialStates().begin()]);
this->samplingModel->instantiate(point);
valueInBoundOfFormula=this->valueIsInBoundOfFormula(this->samplingModel->computeReachabilityProbabilities()[*this->samplingModel->getModel()->getInitialStates().begin()]);
}
if(valueInBoundOfFormula){
@ -827,15 +819,12 @@ namespace storm {
template<typename ParametricType, typename ConstantType>
bool SparseDtmcRegionModelChecker<ParametricType, ConstantType>::checkApproximativeProbabilities(ParameterRegion& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds) {
STORM_LOG_THROW(this->hasOnlyLinearFunctions, storm::exceptions::UnexpectedException, "Tried to perform approximative method (only applicable if all functions are linear), but there are nonlinear functions.");
//build the mdp and a reachability formula and create a modelchecker
std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now();
buildMdpForApproximation(region);
this->approximationModel->instantiate(region);
std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now();
this->timeMDPBuild += timeMDPBuildEnd-timeMDPBuildStart;
std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr);
storm::modelchecker::SparseMdpPrctlModelChecker<ConstantType> modelChecker(*this->approxMdp);
if (region.getCheckResult()==RegionCheckResult::UNKNOWN){
//Sampling to get a hint of whether we should try to prove ALLSAT or ALLVIOLATED
@ -881,13 +870,11 @@ namespace storm {
//However, the second iteration might not be performed if the first one proved ALLSAT or ALLVIOLATED
storm::logic::OptimalityType opType=firstOpType;
for(int i=1; i<=2; ++i){
//perform model checking on the mdp
std::unique_ptr<storm::modelchecker::CheckResult> resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType);
//check if the approximation yielded a conclusive result
if(opType==storm::logic::OptimalityType::Maximize){
upperBounds = resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
if(valueIsInBoundOfFormula(upperBounds[*this->approxMdp->getInitialStates().begin()])){
upperBounds = this->approximationModel->computeReachabilityProbabilities(opType);
if(valueIsInBoundOfFormula(upperBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Less) ||
(this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::LessEqual)){
region.setCheckResult(RegionCheckResult::ALLSAT);
@ -905,8 +892,8 @@ namespace storm {
opType=storm::logic::OptimalityType::Minimize;
}
else if(opType==storm::logic::OptimalityType::Minimize){
lowerBounds = resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
if(valueIsInBoundOfFormula(lowerBounds[*this->approxMdp->getInitialStates().begin()])){
lowerBounds = this->approximationModel->computeReachabilityProbabilities(opType);
if(valueIsInBoundOfFormula(lowerBounds[*this->approximationModel->getModel()->getInitialStates().begin()])){
if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Greater) ||
(this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::GreaterEqual)){
region.setCheckResult(RegionCheckResult::ALLSAT);
@ -924,13 +911,14 @@ namespace storm {
opType=storm::logic::OptimalityType::Maximize;
}
}
//if we reach this point than the result is still inconclusive.
return false;
}
/*
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::buildMdpForApproximation(const ParameterRegion& region) {
//instantiate the substitutions for the given region
std::vector<std::map<VariableType, CoefficientType>> substitutions(this->approxMdpSubstitutions.size());
for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->approxMdpSubstitutions.size(); ++substitutionIndex){
@ -948,133 +936,26 @@ namespace storm {
}
}
std::vector<std::pair<typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, ConstantType>> entryVector;
entryVector.reserve(this->approxMdpMapping.size());
//now put the values into the mdp matrix
for( std::tuple<ParametricType, typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, size_t>& mappingTriple : this->approxMdpMapping){
std::get<1>(mappingTriple).setValue(storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)])
)
);
entryVector.emplace_back(std::get<1>(mappingTriple), storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)])));
}
}
//DELETEME
template<typename ParametricType, typename ConstantType>
storm::models::sparse::Mdp<ConstantType> SparseDtmcRegionModelChecker<ParametricType, ConstantType>::buildMdpForApproximation2(const ParameterRegion& region) {
//We are going to build a (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter
//The matrix (and the Choice labeling)
//the matrix this->sparseTransitions might have empty rows where states have been eliminated.
//The MDP matrix should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly.
//These changes are computed in advance
std::vector<storm::storage::sparse::state_type> newStateIndexMap(this->sparseTransitions.getRowCount(), this->sparseTransitions.getRowCount()); //initialize with some illegal index to easily check if a transition leads to an unselected state
storm::storage::sparse::state_type newStateIndex=0;
for(auto const& state : this->subsystem){
newStateIndexMap[state]=newStateIndex;
++newStateIndex;
std::chrono::high_resolution_clock::time_point timeMDPInsertingStart = std::chrono::high_resolution_clock::now();
for(std::pair<typename storm::storage::MatrixEntry<storm::storage::sparse::state_type, ConstantType>&, ConstantType>& entry : entryVector ){
entry.first.setValue(entry.second);
}
//We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead
storm::storage::sparse::state_type numStates=newStateIndex+2;
storm::storage::sparse::state_type targetState=numStates-2;
storm::storage::sparse::state_type sinkState= numStates-1;
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(0, numStates, 0, true, true, numStates);
//std::vector<boost::container::flat_set<uint_fast64_t>> choiceLabeling;
//fill in the data row by row
storm::storage::sparse::state_type currentRow=0;
for(storm::storage::sparse::state_type oldStateIndex : this->subsystem){
//we first go through the row to find out a) which parameter occur in that row and b) at which point do we have to insert the selfloop
storm::storage::sparse::state_type selfloopIndex=0;
std::set<VariableType> occurringParameters;
for(auto const& entry: this->sparseTransitions.getRow(oldStateIndex)){
storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringParameters);
if(entry.getColumn()<=oldStateIndex){
if(entry.getColumn()==oldStateIndex){
//there already is a selfloop so we do not have to add one.
selfloopIndex=numStates; // --> selfloop will never be inserted
}
else {
++selfloopIndex;
}
}
STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=this->sparseTransitions.getRowCount(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated.");
}
storm::utility::regions::gatherOccurringVariables(this->oneStepProbabilities[oldStateIndex], occurringParameters);
//we now add a row for every combination of lower and upper bound of the variables
auto const& substitutions = region.getVerticesOfRegion(occurringParameters);
STORM_LOG_ASSERT(!substitutions.empty(), "there are no substitutions, i.e., no vertices of the current region. this should not be possible");
matrixBuilder.newRowGroup(currentRow);
for(size_t substitutionIndex=0; substitutionIndex< substitutions.size(); ++substitutionIndex){
ConstantType missingProbability = storm::utility::one<ConstantType>();
if(selfloopIndex==0){ //add the selfloop first.
matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero<ConstantType>());
selfloopIndex=numStates; // --> selfloop will never be inserted again
}
for(auto const& entry : this->sparseTransitions.getRow(oldStateIndex)){
ConstantType value = storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(entry.getValue(),substitutions[substitutionIndex])
);
missingProbability-=value;
storm::storage::sparse::state_type column = newStateIndexMap[entry.getColumn()];
matrixBuilder.addNextValue(currentRow, column, value);
--selfloopIndex;
if(selfloopIndex==0){ //add the selfloop now
matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero<ConstantType>());
selfloopIndex=numStates; // --> selfloop will never be inserted again
}
}
if(!this->parametricTypeComparator.isZero(this->oneStepProbabilities[oldStateIndex])){ //transition to target state
ConstantType value = storm::utility::regions::convertNumber<CoefficientType,ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType,ConstantType>(this->oneStepProbabilities[oldStateIndex],substitutions[substitutionIndex])
);
missingProbability-=value;
matrixBuilder.addNextValue(currentRow, targetState, value);
}
if(!this->constantTypeComparator.isZero(missingProbability)){ //transition to sink state
STORM_LOG_THROW((missingProbability> -storm::utility::regions::convertNumber<double, ConstantType>(storm::settings::generalSettings().getPrecision())), storm::exceptions::UnexpectedException, "The missing probability is negative.");
matrixBuilder.addNextValue(currentRow, sinkState, missingProbability);
}
//boost::container::flat_set<uint_fast64_t> label;
//label.insert(substitutionIndex);
//choiceLabeling.emplace_back(label);
++currentRow;
}
}
//add self loops on the additional states (i.e., target and sink)
//boost::container::flat_set<uint_fast64_t> labelAll;
//labelAll.insert(substitutions.size());
matrixBuilder.newRowGroup(currentRow);
matrixBuilder.addNextValue(currentRow, targetState, storm::utility::one<ConstantType>());
// choiceLabeling.emplace_back(labelAll);
++currentRow;
matrixBuilder.newRowGroup(currentRow);
matrixBuilder.addNextValue(currentRow, sinkState, storm::utility::one<ConstantType>());
// choiceLabeling.emplace_back(labelAll);
++currentRow;
//The labeling
storm::models::sparse::StateLabeling stateLabeling(numStates);
storm::storage::BitVector initLabel(numStates, false);
initLabel.set(newStateIndexMap[this->initialState], true);
stateLabeling.addLabel("init", std::move(initLabel));
storm::storage::BitVector targetLabel(numStates, false);
targetLabel.set(numStates-2, true);
stateLabeling.addLabel("target", std::move(targetLabel));
storm::storage::BitVector sinkLabel(numStates, false);
sinkLabel.set(numStates-1, true);
stateLabeling.addLabel("sink", std::move(sinkLabel));
// The MDP
boost::optional<std::vector<ConstantType>> noStateRewards;
boost::optional<storm::storage::SparseMatrix<ConstantType>> noTransitionRewards;
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;
return storm::models::sparse::Mdp<ConstantType>(matrixBuilder.build(), std::move(stateLabeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
std::chrono::high_resolution_clock::time_point timeMDPInsertingEnd = std::chrono::high_resolution_clock::now();
this->timeMDPInserting += timeMDPInsertingEnd-timeMDPInsertingStart;
}
*/
@ -1214,19 +1095,12 @@ namespace storm {
std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ...
std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(timeOverall);
std::size_t subsystemTransitions = this->sparseTransitions.getNonzeroEntryCount();
for(auto const& transition : this->oneStepProbabilities){
if(!this->parametricTypeComparator.isZero(transition)){
++subsystemTransitions;
}
}
uint_fast64_t numOfSolvedRegions= this->numOfRegionsExistsBoth + this->numOfRegionsAllSat + this->numOfRegionsAllViolated;
outstream << std::endl << "Statistics Region Model Checker Statistics:" << std::endl;
outstream << "-----------------------------------------------" << std::endl;
outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl;
outstream << "Reduced model: " << this->subsystem.getNumberOfSetBits() << " states, " << subsystemTransitions << " transitions" << std::endl;
outstream << "Simplified model: " << this->simplifiedModel->getNumberOfStates() << " states, " << this->simplifiedModel->getNumberOfTransitions() << " transitions" << std::endl;
outstream << "Formula: " << *this->probabilityOperatorFormula << std::endl;
outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl;
outstream << "Number of checked regions: " << this->numOfCheckedRegions << std::endl;

69
src/modelchecker/reachability/SparseDtmcRegionModelChecker.h

@ -9,7 +9,9 @@
#include "src/utility/constants.h"
#include "utility/regions.h"
#include "src/solver/Smt2SmtSolver.h"
#include "SparseDtmcEliminationModelChecker.h"
#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
namespace storm {
namespace modelchecker {
@ -18,6 +20,8 @@ namespace storm {
class SparseDtmcRegionModelChecker {
public:
//The type of variables and bounds depends on the template arguments
typedef typename std::conditional<(std::is_same<ParametricType,storm::RationalFunction>::value), storm::Variable,std::nullptr_t>::type VariableType;
typedef typename std::conditional<(std::is_same<ParametricType,storm::RationalFunction>::value), storm::Coefficient,std::nullptr_t>::type CoefficientType;
@ -34,6 +38,9 @@ namespace storm {
ALLVIOLATED /*!< the formula is violated for all parameters in the given region */
};
class ApproximationModel;
class SamplingModel;
class ParameterRegion{
public:
@ -202,21 +209,13 @@ namespace storm {
* eliminates all states for which the outgoing transitions are constant.
* Also checks whether the non constant functions are linear
*/
void eliminateStatesConstSucc(
storm::storage::BitVector& subsys,
FlexibleMatrix& flexTransitions,
FlexibleMatrix& flexBackwardTransitions,
std::vector<ParametricType>& oneStepProbs,
bool& allFunctionsAreLinear,
storm::storage::sparse::state_type const& initState
);
void computeSimplifiedModel(storm::storage::BitVector const& targetStates);
enum class TypeOfBound {
LOWER,
UPPER
};
/*!
* Initializes a Sample-Model that can be used to get the probability result for a certain parameter evaluation.
* Initializes an Approximation-Model that can be used to approximate the reachability probabilities.
*/
void initializeSampleAndApproxModel();
/*!
* Initializes a DTMC that can be used to get the probability result for a certain parameter evaluation.
@ -236,8 +235,7 @@ namespace storm {
* @param flexTransitions the transitions of the pDTMC
* @param oneStepProbs the probabilities to move to a target state
* @param initState the initial state of the pDtmc
*/
void initializeSampleDtmcAndApproxMdp(
void initializeSampleDtmcAndApproxMdp2(
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> & sampleDtmc,
std::vector<std::pair<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&>> & sampleDtmcMapping,
std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> & approxMdp,
@ -248,19 +246,10 @@ namespace storm {
std::vector<ParametricType> const& oneStepProbs,
storm::storage::sparse::state_type const& initState
);
*/
ParametricType getReachProbFunction();
//Computes the reachability probability function by performing state elimination
ParametricType computeReachProbFunction(
storm::storage::BitVector const& subsys,
FlexibleMatrix const& flexTransitions,
FlexibleMatrix const& flexBackwardTransitions,
storm::storage::SparseMatrix<ParametricType> const& spTransitions,
storm::storage::SparseMatrix<ParametricType> const& spBackwardTransitions,
std::vector<ParametricType> const& oneStepProbs,
storm::storage::sparse::state_type const& initState
);
//initializes the given solver which can later be used to give an exact result regarding the whole model.
void initializeSMTSolver(std::shared_ptr<storm::solver::Smt2SmtSolver>& solver, ParametricType const& reachProbFunction, storm::logic::ProbabilityOperatorFormula const& formula);
@ -339,21 +328,18 @@ namespace storm {
//the currently specified formula
std::unique_ptr<storm::logic::ProbabilityOperatorFormula> probabilityOperatorFormula;
// The ingredients of the model where constant transitions have been eliminated as much as possible
// the probability matrix
FlexibleMatrix flexibleTransitions;
storm::storage::SparseMatrix<ParametricType> sparseTransitions;
//the corresponding backward transitions
FlexibleMatrix flexibleBackwardTransitions;
storm::storage::SparseMatrix<ParametricType> sparseBackwardTransitions;
// the propabilities to go to a state with probability 1 in one step (belongs to flexibleTransitions)
std::vector<ParametricType> oneStepProbabilities;
// the initial state
storm::storage::sparse::state_type initialState;
// the set of states that have not been eliminated
storm::storage::BitVector subsystem;
// the original model after states with constant transitions have been eliminated
std::shared_ptr<storm::models::sparse::Dtmc<ParametricType>> simplifiedModel;
// a flag that is true if there are only linear functions at transitions of the model
bool hasOnlyLinearFunctions;
// the model that can be instantiated to check the value at a certain point
std::shared_ptr<SamplingModel> samplingModel;
// the model that is used to approximate the probability values
std::shared_ptr<ApproximationModel> approximationModel;
/*
// the dtmc that can be instantiated to check the value at a certain point
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> sampleDtmc;
// a vector that links entries of the dtmc matrix with the corresponding transition functions (for fast instantiation)
@ -364,10 +350,11 @@ namespace storm {
std::vector<std::tuple<ParametricType, storm::storage::MatrixEntry<storm::storage::sparse::state_type,ConstantType>&, std::size_t>> approxMdpMapping;
// stores the different substitutions of the variables
std::vector<std::map<VariableType, TypeOfBound>> approxMdpSubstitutions;
*/
// The function for the reachability probability in the initial state
ParametricType reachProbFunction;
bool isReachProbFunctionComputed;
bool isResultConstant;
// runtimes and other information for statistics.

170
src/modelchecker/region/ApproximationModel.cpp

@ -0,0 +1,170 @@
/*
* File: ApproximationModel.cpp
* Author: tim
*
* Created on August 7, 2015, 9:29 AM
*/
#include "src/modelchecker/region/ApproximationModel.h"
#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h"
#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "exceptions/UnexpectedException.h"
namespace storm {
namespace modelchecker {
template<typename ParametricType, typename ConstantType>
SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::ApproximationModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel) : mapping(), evaluationTable(), substitutions() {
// Run through the rows of the original model and obtain
// (1) the different substitutions (this->substitutions) and the substitution used for every row
std::vector<std::size_t> rowSubstitutions;
// (2) the set of distinct pairs <Function, Substitution>
std::set<std::pair<ParametricType, std::size_t>> distinctFuncSubs;
// (3) the MDP matrix with some dummy entries
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates());
typename storm::storage::SparseMatrix<ConstantType>::index_type approxModelRow=0;
uint_fast64_t numOfNonConstEntries=0;
for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
matrixBuilder.newRowGroup(approxModelRow);
// For (1): 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(row)){
storm::utility::regions::gatherOccurringVariables(entry.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)
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
//(00...0 = lower bounds for all parameters, 11...1 = upper bounds for all parameters)
std::map<VariableType, TypeOfBound> 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));
}
++parameterIndex;
}
//Find the current substitution in this->substitutions (add it if we see this substitution the first time)
std::size_t substitutionIndex = std::find(this->substitutions.begin(),this->substitutions.end(), currSubstitution) - this->substitutions.begin();
if(substitutionIndex==this->substitutions.size()){
this->substitutions.emplace_back(std::move(currSubstitution));
}
rowSubstitutions.push_back(substitutionIndex);
//Run again through the row and...
//For (2): add pair <function, substitution> for the occuring functions and the current substitution
//For (3): add a row with dummy entries for the current substitution
//Note that this is still executed even if no variables occur. However, we will not put any <function, substitution> pair into distninctFuncSubs if substitution is empty.
ConstantType dummyEntry=storm::utility::one<ConstantType>();
for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
if(!this->parametricTypeComparator.isConstant(entry.getValue())){
distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex));
++numOfNonConstEntries;
}
matrixBuilder.addNextValue(approxModelRow, entry.getColumn(), dummyEntry);
dummyEntry=storm::utility::zero<ConstantType>();
}
++approxModelRow;
}
}
//Obtain other model ingredients and the approximation model itself
storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
boost::optional<std::vector<ConstantType>> noStateRewards;
boost::optional<storm::storage::SparseMatrix<ConstantType>> noTransitionRewards;
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;
this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
//Get the evaluation table. Note that it remains sorted due to the fact that distinctFuncSubs is sorted
this->evaluationTable.reserve(distinctFuncSubs.size());
for(auto const& funcSub : distinctFuncSubs){
this->evaluationTable.emplace_back(funcSub.first, funcSub.second, storm::utility::zero<ConstantType>());
}
//Fill in the entries for the mapping
this->mapping.reserve(numOfNonConstEntries);
approxModelRow=0;
for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
for (; approxModelRow<this->model->getTransitionMatrix().getRowGroupIndices()[row+1]; ++approxModelRow){
auto appEntry = this->model->getTransitionMatrix().getRow(approxModelRow).begin();
for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){
if(this->parametricTypeComparator.isConstant(parEntry.getValue())){
appEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart<ParametricType, ConstantType>(parEntry.getValue())));
}
else {
std::tuple<ParametricType, std::size_t, ConstantType> searchedTuple(parEntry.getValue(), rowSubstitutions[approxModelRow], storm::utility::zero<ConstantType>());
auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedTuple);
STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the evaluationTable. Either the table is missing that tuple or it is not sorted properly");
mapping.emplace_back(std::make_pair(&(std::get<2>(*tableIt)), appEntry));
}
++appEntry;
}
}
}
}
template<typename ParametricType, typename ConstantType>
SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::~ApproximationModel() {
//Intentionally left empty
}
template<typename ParametricType, typename ConstantType>
std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::getModel() const {
return this->model;
}
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::instantiate(const ParameterRegion& region) {
//Instantiate the substitutions
std::vector<std::map<VariableType, CoefficientType>> instantiatedSubs(this->substitutions.size());
for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->substitutions.size(); ++substitutionIndex){
for(std::pair<VariableType, TypeOfBound> const& sub : this->substitutions[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->evaluationTable){
std::get<2>(tableEntry)=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(
std::get<0>(tableEntry),
instantiatedSubs[std::get<1>(tableEntry)]
)
);
}
//write the instantiated values to the matrix according to the mapping
for(auto& mappingPair : this->mapping){
mappingPair.second->setValue(*(mappingPair.first));
}
}
template<typename ParametricType, typename ConstantType>
std::vector<ConstantType> const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::computeReachabilityProbabilities(storm::logic::OptimalityType const& optimalityType) {
std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr);
storm::modelchecker::SparseMdpPrctlModelChecker<ConstantType> modelChecker(*this->model);
//perform model checking on the mdp
std::unique_ptr<storm::modelchecker::CheckResult> resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula, false, optimalityType);
return resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
}
#ifdef STORM_HAVE_CARL
template class SparseDtmcRegionModelChecker<storm::RationalFunction, double>::ApproximationModel;
#endif
}
}

82
src/modelchecker/region/ApproximationModel.h

@ -0,0 +1,82 @@
/*
* File: ApproximationModel.h
* Author: tim
*
* Created on August 7, 2015, 9:29 AM
*/
#ifndef STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H
#define STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H
#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
namespace storm {
namespace modelchecker {
template<typename ParametricType, typename ConstantType>
class SparseDtmcRegionModelChecker;
template<typename ParametricType, typename ConstantType>
class SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel{
public:
typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType VariableType;
typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType CoefficientType;
ApproximationModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel);
virtual ~ApproximationModel();
/*!
* returns the underlying model
*/
std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> const& getModel() const;
/*!
* Instantiates the underlying model according to the given region
*/
void instantiate(ParameterRegion const& region);
/*!
* Returns the approximated reachability probabilities for every state.
* Undefined behavior if model has not been instantiated first!
* @param optimalityType Use MAXIMIZE to get upper bounds or MINIMIZE to get lower bounds
*/
std::vector<ConstantType> const& computeReachabilityProbabilities(storm::logic::OptimalityType const& optimalityType);
private:
enum class TypeOfBound {
LOWER,
UPPER
};
//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>> mapping;
//Vector has one entry for
//(every distinct, non-constant function that occurs somewhere in the model) x (the required combinations of lower and upper bounds of the region)
//The second entry represents a substitution as an index in the substitutions vector
//The third entry should contain the result when evaluating the function in the first entry, regarding the substitution given by the second entry.
std::vector<std::tuple<ParametricType, std::size_t, ConstantType>> evaluationTable;
//Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region)
std::vector<std::map<VariableType, TypeOfBound>> substitutions;
//The Model with which we work
std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> model;
// comparators that can be used to compare constants.
storm::utility::ConstantsComparator<ParametricType> parametricTypeComparator;
storm::utility::ConstantsComparator<ConstantType> constantTypeComparator;
};
}
}
#endif /* STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H */

107
src/modelchecker/region/SamplingModel.cpp

@ -0,0 +1,107 @@
/*
* File: SamplingModel.cpp
* Author: tim
*
* Created on August 7, 2015, 9:31 AM
*/
#include "src/modelchecker/region/SamplingModel.h"
#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "exceptions/UnexpectedException.h"
namespace storm {
namespace modelchecker {
template<typename ParametricType, typename ConstantType>
SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::SamplingModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel) : mapping(), evaluationTable(){
// Run through the rows of the original model and obtain the set of distinct functions as well as a matrix with dummy entries
std::set<ParametricType> functionSet;
storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount());
uint_fast64_t numOfNonConstEntries=0;
for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
ConstantType dummyEntry=storm::utility::one<ConstantType>();
for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
if(!this->parametricTypeComparator.isConstant(entry.getValue())){
functionSet.insert(entry.getValue());
++numOfNonConstEntries;
}
matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry);
dummyEntry=storm::utility::zero<ConstantType>();
}
}
//Obtain other model ingredients and the approximation model itself
storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
boost::optional<std::vector<ConstantType>> noStateRewards;
boost::optional<storm::storage::SparseMatrix<ConstantType>> noTransitionRewards;
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;
this->model=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
//Get the evaluation table. Note that it remains sorted due to the fact that functionSet is sorted
this->evaluationTable.reserve(functionSet.size());
for(auto const& func : functionSet){
this->evaluationTable.emplace_back(func, storm::utility::zero<ConstantType>());
}
//Fill in the entries for the mapping
this->mapping.reserve(numOfNonConstEntries);
auto samEntry = this->model->getTransitionMatrix().begin();
for(auto const& parEntry : parametricModel.getTransitionMatrix()){
if(this->parametricTypeComparator.isConstant(parEntry.getValue())){
samEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart<ParametricType, ConstantType>(parEntry.getValue())));
}
else {
std::pair<ParametricType, ConstantType> searchedPair(parEntry.getValue(), storm::utility::zero<ConstantType>());
auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair);
STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly");
mapping.emplace_back(std::make_pair(&(tableIt->second), samEntry));
}
++samEntry;
}
}
template<typename ParametricType, typename ConstantType>
SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::~SamplingModel() {
//Intentionally left empty
}
template<typename ParametricType, typename ConstantType>
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::getModel() const {
return this->model;
}
template<typename ParametricType, typename ConstantType>
void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::instantiate(std::map<VariableType, CoefficientType>const& point) {
//write entries into evaluation table
for(auto& tableEntry : this->evaluationTable){
tableEntry.second=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
storm::utility::regions::evaluateFunction<ParametricType, ConstantType>(
tableEntry.first,
point
)
);
}
//write the instantiated values to the matrix according to the mapping
for(auto& mappingPair : this->mapping){
mappingPair.second->setValue(*(mappingPair.first));
}
}
template<typename ParametricType, typename ConstantType>
std::vector<ConstantType> const& SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::computeReachabilityProbabilities() {
std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr);
storm::modelchecker::SparseDtmcPrctlModelChecker<ConstantType> modelChecker(*this->model);
//perform model checking on the mdp
std::unique_ptr<storm::modelchecker::CheckResult> resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula);
return resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
}
#ifdef STORM_HAVE_CARL
template class SparseDtmcRegionModelChecker<storm::RationalFunction, double>::SamplingModel;
#endif
}
}

72
src/modelchecker/region/SamplingModel.h

@ -0,0 +1,72 @@
/*
* File: SamplingModel.h
* Author: tim
*
* Created on August 7, 2015, 9:31 AM
*/
#ifndef STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H
#define STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H
#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h"
namespace storm {
namespace modelchecker{
template<typename ParametricType, typename ConstantType>
class SparseDtmcRegionModelChecker;
template<typename ParametricType, typename ConstantType>
class SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel {
public:
typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::VariableType VariableType;
typedef typename SparseDtmcRegionModelChecker<ParametricType, ConstantType>::CoefficientType CoefficientType;
SamplingModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel);
virtual ~SamplingModel();
/*!
* returns the underlying model
*/
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> const& getModel() const;
/*!
* Instantiates the underlying model according to the given point
*/
void instantiate(std::map<VariableType, CoefficientType>const& point);
/*!
* Returns the reachability probabilities for every state according to the current instantiation.
* Undefined behavior if model has not been instantiated first!
*
* @param optimalityType Use MAXIMIZE to get upper bounds or MINIMIZE to get lower bounds
*/
std::vector<ConstantType> const& computeReachabilityProbabilities();
private:
//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>> mapping;
//Vector has one entry for every distinct, non-constant function that occurs somewhere in the model.
//The second entry should contain the result when evaluating the function in the first entry.
std::vector<std::pair<ParametricType, ConstantType>> evaluationTable;
//The model with which we work
std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> model;
// comparators that can be used to compare constants.
storm::utility::ConstantsComparator<ParametricType> parametricTypeComparator;
storm::utility::ConstantsComparator<ConstantType> constantTypeComparator;
};
}
}
#endif /* STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H */

5
src/utility/regions.cpp

@ -184,6 +184,11 @@ namespace storm {
return function.evaluate(point);
}
template<>
typename storm::modelchecker::SparseDtmcRegionModelChecker<storm::RationalFunction,double>::CoefficientType getConstantPart<storm::RationalFunction, double>(storm::RationalFunction const& function){
return function.constantPart();
}
template<>
bool functionIsLinear<storm::RationalFunction>(storm::RationalFunction const& function){
// Note: At this moment there is no function in carl for rationalFunctions.

10
src/utility/regions.h

@ -114,9 +114,19 @@ namespace storm {
template<typename VariableType>
std::string getVariableName(VariableType variable);
/*
* evaluates the given function at the given point and returns the result
*/
template<typename ParametricType, typename ConstantType>
typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType evaluateFunction(ParametricType const& function, std::map<typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::VariableType, typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType> const& point);
/*!
* retrieves the constant part of the given function.
* If the function is constant, then the result is the same value as the original function
*/
template<typename ParametricType, typename ConstantType>
typename storm::modelchecker::SparseDtmcRegionModelChecker<ParametricType,ConstantType>::CoefficientType getConstantPart(ParametricType const& function);
/*!
* Returns true if the function is rational. Note that the function might be simplified.
*/

Loading…
Cancel
Save