From 28326b14e41ee837d2b7e7d1300e3a31fd5a3155 Mon Sep 17 00:00:00 2001 From: TimQu Date: Sun, 9 Aug 2015 19:10:27 +0200 Subject: [PATCH] more efficient instantiation of matrices. Former-commit-id: 0f4c36f2f21db4455597858b925e190a15f89eb2 --- .../SparseDtmcRegionModelChecker.cpp | 458 +++++++----------- .../SparseDtmcRegionModelChecker.h | 69 ++- .../region/ApproximationModel.cpp | 170 +++++++ src/modelchecker/region/ApproximationModel.h | 82 ++++ src/modelchecker/region/SamplingModel.cpp | 107 ++++ src/modelchecker/region/SamplingModel.h | 72 +++ src/utility/regions.cpp | 5 + src/utility/regions.h | 10 + 8 files changed, 640 insertions(+), 333 deletions(-) create mode 100644 src/modelchecker/region/ApproximationModel.cpp create mode 100644 src/modelchecker/region/ApproximationModel.h create mode 100644 src/modelchecker/region/SamplingModel.cpp create mode 100644 src/modelchecker/region/SamplingModel.h diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp index 9619f4c6d..66378024b 100644 --- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp +++ b/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(new storm::logic::ProbabilityOperatorFormula(formula.asStateFormula().asProbabilityOperatorFormula())); storm::logic::EventuallyFormula const& eventuallyFormula = this->probabilityOperatorFormula->getSubformula().asPathFormula().asEventuallyFormula(); std::unique_ptr 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 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() : storm::utility::one(); - 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 submatrix = model.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates); - storm::storage::SparseMatrix 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 - void SparseDtmcRegionModelChecker::eliminateStatesConstSucc( - storm::storage::BitVector& subsys, - FlexibleMatrix& flexTransitions, - FlexibleMatrix& flexBackwardTransitions, - std::vector& oneStepProbs, - bool& allFunctionsAreLinear, - storm::storage::sparse::state_type const& initialState) { + void SparseDtmcRegionModelChecker::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 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() : storm::utility::one(); + 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 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 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> 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(entry.getValue()); + this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear(entry.getValue()); } } - if(!this->parametricTypeComparator.isConstant(oneStepProbs[state])){ + if(!this->parametricTypeComparator.isConstant(oneStepProbabilities[state])){ onlyConstantOutgoingTransitions=false; - allFunctionsAreLinear &= storm::utility::regions::functionIsLinear(oneStepProbs[state]); + this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear(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 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 matrixBuilder(numStates, numStates); + for(storm::storage::sparse::state_type oldStateIndex : subsystem){ + ParametricType valueToSinkState=storm::utility::regions::getNewFunction(storm::utility::one()); + //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()); + matrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one()); + //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> noStateRewards; + boost::optional> noTransitionRewards; + boost::optional>> noChoiceLabeling; + // the final model + this->simplifiedModel = std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); } template - void SparseDtmcRegionModelChecker::initializeSampleDtmcAndApproxMdp( + void SparseDtmcRegionModelChecker::initializeSampleAndApproxModel(){ + + + + //now create the models + this->approximationModel=std::make_shared(*this->simplifiedModel); + this->samplingModel=std::make_shared(*this->simplifiedModel); + + } + /* + template + void SparseDtmcRegionModelChecker::initializeSampleDtmcAndApproxMdp2( std::shared_ptr>& sampleDtmc, std::vector&>> & sampleDtmcMapping, std::shared_ptr> & approxMdp, @@ -344,7 +391,6 @@ namespace storm { std::vector const& oneStepProbs, storm::storage::sparse::state_type const& initState) { - std::set 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>(dtmcMatrixBuilder.build(), std::move(dtmcStateLabeling), std::move(dtmcNoStateRewards), std::move(dtmcNoTransitionRewards), std::move(dtmcNoChoiceLabeling)); approxMdp=std::make_shared>(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&>>(); sampleDtmcMapping.reserve(sampleDtmc->getTransitionMatrix().getEntryCount()); @@ -532,12 +574,29 @@ namespace storm { } } } - + */ template ParametricType SparseDtmcRegionModelChecker::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 backwardTransitions(this->simplifiedModel->getBackwardTransitions()); + std::vector oneStepProbabilities(this->simplifiedModel->getNumberOfStates()-2, storm::utility::zero()); + for(auto const& entry : backwardTransitions.getRow(*(this->simplifiedModel->getStates("target").begin()))){ + if(entry.getColumn()simplifiedModel->getStates("target") | this->simplifiedModel->getStates("sink")); + backwardTransitions=backwardTransitions.getSubmatrix(false,maybeStates,maybeStates); + storm::storage::SparseMatrix 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> noStateRewards; + std::vector 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 - ParametricType SparseDtmcRegionModelChecker::computeReachProbFunction( - storm::storage::BitVector const& subsys, - FlexibleMatrix const& flexTransitions, - FlexibleMatrix const& flexBackwardTransitions, - storm::storage::SparseMatrix const& spTransitions, - storm::storage::SparseMatrix const& spBackwardTransitions, - std::vector 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 workingCopyOneStepProbs(oneStepProbs); - - storm::storage::BitVector initialStates(subsys.size(),false); - initialStates.set(initState, true); - std::vector statePriorities = this->eliminationModelChecker.getStatePriorities(spTransitions, spBackwardTransitions, initialStates, workingCopyOneStepProbs); - boost::optional> 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 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 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 void SparseDtmcRegionModelChecker::initializeSMTSolver(std::shared_ptr& solver, const ParametricType& reachProbFunc, const storm::logic::ProbabilityOperatorFormula& formula) { @@ -788,18 +790,8 @@ namespace storm { valueInBoundOfFormula = this->valueIsInBoundOfFormula(storm::utility::regions::evaluateFunction(getReachProbFunction(), point)); } else{ - //put the values into the dtmc matrix - for( std::pair&>& mappingPair : this->sampleDtmcMapping){ - mappingPair.second.setValue(storm::utility::regions::convertNumber( - storm::utility::regions::evaluateFunction(mappingPair.first, point) - ) - ); - } - std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); - storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); - storm::modelchecker::SparseDtmcPrctlModelChecker modelChecker(*this->sampleDtmc); - std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false); - valueInBoundOfFormula=this->valueIsInBoundOfFormula(resultPtr->asExplicitQuantitativeCheckResult().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 bool SparseDtmcRegionModelChecker::checkApproximativeProbabilities(ParameterRegion& region, std::vector& lowerBounds, std::vector& 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 targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); - storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); - storm::modelchecker::SparseMdpPrctlModelChecker 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 resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType); //check if the approximation yielded a conclusive result if(opType==storm::logic::OptimalityType::Maximize){ - upperBounds = resultPtr->asExplicitQuantitativeCheckResult().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().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 void SparseDtmcRegionModelChecker::buildMdpForApproximation(const ParameterRegion& region) { + //instantiate the substitutions for the given region std::vector> substitutions(this->approxMdpSubstitutions.size()); for(uint_fast64_t substitutionIndex=0; substitutionIndexapproxMdpSubstitutions.size(); ++substitutionIndex){ @@ -948,133 +936,26 @@ namespace storm { } } + std::vector&, ConstantType>> entryVector; + entryVector.reserve(this->approxMdpMapping.size()); //now put the values into the mdp matrix for( std::tuple&, size_t>& mappingTriple : this->approxMdpMapping){ - std::get<1>(mappingTriple).setValue(storm::utility::regions::convertNumber( - storm::utility::regions::evaluateFunction(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)]) - ) - ); + entryVector.emplace_back(std::get<1>(mappingTriple), storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)]))); } - } - - //DELETEME - template - storm::models::sparse::Mdp SparseDtmcRegionModelChecker::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 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&, 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 matrixBuilder(0, numStates, 0, true, true, numStates); - //std::vector> 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 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(); - if(selfloopIndex==0){ //add the selfloop first. - matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero()); - selfloopIndex=numStates; // --> selfloop will never be inserted again - } - for(auto const& entry : this->sparseTransitions.getRow(oldStateIndex)){ - ConstantType value = storm::utility::regions::convertNumber( - storm::utility::regions::evaluateFunction(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()); - 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( - storm::utility::regions::evaluateFunction(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(storm::settings::generalSettings().getPrecision())), storm::exceptions::UnexpectedException, "The missing probability is negative."); - matrixBuilder.addNextValue(currentRow, sinkState, missingProbability); - } - //boost::container::flat_set 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 labelAll; - //labelAll.insert(substitutions.size()); - matrixBuilder.newRowGroup(currentRow); - matrixBuilder.addNextValue(currentRow, targetState, storm::utility::one()); - // choiceLabeling.emplace_back(labelAll); - ++currentRow; - matrixBuilder.newRowGroup(currentRow); - matrixBuilder.addNextValue(currentRow, sinkState, storm::utility::one()); - // 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> noStateRewards; - boost::optional> noTransitionRewards; - boost::optional>> noChoiceLabeling; - - return storm::models::sparse::Mdp(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(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; diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h index f2ed8f521..7f58754f3 100644 --- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h +++ b/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::value), storm::Variable,std::nullptr_t>::type VariableType; typedef typename std::conditional<(std::is_same::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& 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> & sampleDtmc, std::vector&>> & sampleDtmcMapping, std::shared_ptr> & approxMdp, @@ -248,19 +246,10 @@ namespace storm { std::vector 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 const& spTransitions, - storm::storage::SparseMatrix const& spBackwardTransitions, - std::vector 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& solver, ParametricType const& reachProbFunction, storm::logic::ProbabilityOperatorFormula const& formula); @@ -339,21 +328,18 @@ namespace storm { //the currently specified formula std::unique_ptr probabilityOperatorFormula; - // The ingredients of the model where constant transitions have been eliminated as much as possible - // the probability matrix - FlexibleMatrix flexibleTransitions; - storm::storage::SparseMatrix sparseTransitions; - //the corresponding backward transitions - FlexibleMatrix flexibleBackwardTransitions; - storm::storage::SparseMatrix sparseBackwardTransitions; - // the propabilities to go to a state with probability 1 in one step (belongs to flexibleTransitions) - std::vector 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> 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; + // the model that is used to approximate the probability values + std::shared_ptr approximationModel; + + /* // the dtmc that can be instantiated to check the value at a certain point std::shared_ptr> 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::size_t>> approxMdpMapping; // stores the different substitutions of the variables std::vector> approxMdpSubstitutions; - + */ // The function for the reachability probability in the initial state ParametricType reachProbFunction; bool isReachProbFunctionComputed; + bool isResultConstant; // runtimes and other information for statistics. diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp new file mode 100644 index 000000000..3e157e62b --- /dev/null +++ b/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 + SparseDtmcRegionModelChecker::ApproximationModel::ApproximationModel(storm::models::sparse::Dtmc 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 rowSubstitutions; + // (2) the set of distinct pairs + std::set> distinctFuncSubs; + // (3) the MDP matrix with some dummy entries + storm::storage::SparseMatrixBuilder matrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates()); + typename storm::storage::SparseMatrix::index_type approxModelRow=0; + uint_fast64_t numOfNonConstEntries=0; + for(typename storm::storage::SparseMatrix::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 occurringVariables; + for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ + storm::utility::regions::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)); + } + ++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 for the occuring functions and the current substitution + //For (3): add a row with dummy entries for the current substitution + //Note that this is still executed even if no variables occur. However, we will not put any pair into distninctFuncSubs if substitution is empty. + ConstantType dummyEntry=storm::utility::one(); + for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ + if(!this->parametricTypeComparator.isConstant(entry.getValue())){ + distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex)); + ++numOfNonConstEntries; + } + matrixBuilder.addNextValue(approxModelRow, entry.getColumn(), dummyEntry); + dummyEntry=storm::utility::zero(); + } + ++approxModelRow; + } + } + + //Obtain other model ingredients and the approximation model itself + storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); + boost::optional> noStateRewards; + boost::optional> noTransitionRewards; + boost::optional>> noChoiceLabeling; + this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + + //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()); + } + + //Fill in the entries for the mapping + this->mapping.reserve(numOfNonConstEntries); + approxModelRow=0; + for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ + for (; approxModelRowmodel->getTransitionMatrix().getRowGroupIndices()[row+1]; ++approxModelRow){ + auto appEntry = this->model->getTransitionMatrix().getRow(approxModelRow).begin(); + for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){ + if(this->parametricTypeComparator.isConstant(parEntry.getValue())){ + appEntry->setValue(storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parEntry.getValue()))); + } + else { + std::tuple searchedTuple(parEntry.getValue(), rowSubstitutions[approxModelRow], storm::utility::zero()); + auto const tableIt= std::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 + SparseDtmcRegionModelChecker::ApproximationModel::~ApproximationModel() { + //Intentionally left empty + } + + template + std::shared_ptr> const& SparseDtmcRegionModelChecker::ApproximationModel::getModel() const { + return this->model; + } + + template + void SparseDtmcRegionModelChecker::ApproximationModel::instantiate(const ParameterRegion& region) { + //Instantiate the substitutions + std::vector> instantiatedSubs(this->substitutions.size()); + for(uint_fast64_t substitutionIndex=0; substitutionIndexsubstitutions.size(); ++substitutionIndex){ + for(std::pair const& sub : this->substitutions[substitutionIndex]){ + 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( + storm::utility::regions::evaluateFunction( + 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 + std::vector const& SparseDtmcRegionModelChecker::ApproximationModel::computeReachabilityProbabilities(storm::logic::OptimalityType const& optimalityType) { + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); + storm::modelchecker::SparseMdpPrctlModelChecker modelChecker(*this->model); + + //perform model checking on the mdp + std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula, false, optimalityType); + return resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); + } + + + +#ifdef STORM_HAVE_CARL + template class SparseDtmcRegionModelChecker::ApproximationModel; +#endif + } +} \ No newline at end of file diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h new file mode 100644 index 000000000..b2da957b9 --- /dev/null +++ b/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 + class SparseDtmcRegionModelChecker; + + template + class SparseDtmcRegionModelChecker::ApproximationModel{ + + public: + + typedef typename SparseDtmcRegionModelChecker::VariableType VariableType; + typedef typename SparseDtmcRegionModelChecker::CoefficientType CoefficientType; + + + ApproximationModel(storm::models::sparse::Dtmc const& parametricModel); + virtual ~ApproximationModel(); + + /*! + * returns the underlying model + */ + std::shared_ptr> 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 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::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> evaluationTable; + + //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region) + std::vector> substitutions; + + //The Model with which we work + std::shared_ptr> model; + + // comparators that can be used to compare constants. + storm::utility::ConstantsComparator parametricTypeComparator; + storm::utility::ConstantsComparator constantTypeComparator; + + }; + } +} + + +#endif /* STORM_MODELCHECKER_REGION_APPROXIMATIONMODEL_H */ + diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp new file mode 100644 index 000000000..7ed31bd37 --- /dev/null +++ b/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 + SparseDtmcRegionModelChecker::SamplingModel::SamplingModel(storm::models::sparse::Dtmc 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 functionSet; + storm::storage::SparseMatrixBuilder matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount()); + uint_fast64_t numOfNonConstEntries=0; + for(typename storm::storage::SparseMatrix::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){ + ConstantType dummyEntry=storm::utility::one(); + for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){ + if(!this->parametricTypeComparator.isConstant(entry.getValue())){ + functionSet.insert(entry.getValue()); + ++numOfNonConstEntries; + } + matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry); + dummyEntry=storm::utility::zero(); + } + } + + //Obtain other model ingredients and the approximation model itself + storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); + boost::optional> noStateRewards; + boost::optional> noTransitionRewards; + boost::optional>> noChoiceLabeling; + this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + + //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()); + } + + //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(storm::utility::regions::getConstantPart(parEntry.getValue()))); + } + else { + std::pair searchedPair(parEntry.getValue(), storm::utility::zero()); + auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair); + STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly"); + mapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); + } + ++samEntry; + } + } + + template + SparseDtmcRegionModelChecker::SamplingModel::~SamplingModel() { + //Intentionally left empty + } + + template + std::shared_ptr> const& SparseDtmcRegionModelChecker::SamplingModel::getModel() const { + return this->model; + } + + template + void SparseDtmcRegionModelChecker::SamplingModel::instantiate(std::mapconst& point) { + //write entries into evaluation table + for(auto& tableEntry : this->evaluationTable){ + tableEntry.second=storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction( + 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 + std::vector const& SparseDtmcRegionModelChecker::SamplingModel::computeReachabilityProbabilities() { + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); + storm::modelchecker::SparseDtmcPrctlModelChecker modelChecker(*this->model); + + //perform model checking on the mdp + std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula); + return resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); + } + +#ifdef STORM_HAVE_CARL + template class SparseDtmcRegionModelChecker::SamplingModel; +#endif + } +} \ No newline at end of file diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h new file mode 100644 index 000000000..ffa50f221 --- /dev/null +++ b/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 + class SparseDtmcRegionModelChecker; + + template + class SparseDtmcRegionModelChecker::SamplingModel { + + public: + + typedef typename SparseDtmcRegionModelChecker::VariableType VariableType; + typedef typename SparseDtmcRegionModelChecker::CoefficientType CoefficientType; + + SamplingModel(storm::models::sparse::Dtmc const& parametricModel); + virtual ~SamplingModel(); + + /*! + * returns the underlying model + */ + std::shared_ptr> const& getModel() const; + + /*! + * Instantiates the underlying model according to the given point + */ + void instantiate(std::mapconst& 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 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::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> evaluationTable; + + //The model with which we work + std::shared_ptr> model; + + // comparators that can be used to compare constants. + storm::utility::ConstantsComparator parametricTypeComparator; + storm::utility::ConstantsComparator constantTypeComparator; + }; + + } +} +#endif /* STORM_MODELCHECKER_REGION_SAMPLINGMODEL_H */ + diff --git a/src/utility/regions.cpp b/src/utility/regions.cpp index 14d261de9..7fffea186 100644 --- a/src/utility/regions.cpp +++ b/src/utility/regions.cpp @@ -184,6 +184,11 @@ namespace storm { return function.evaluate(point); } + template<> + typename storm::modelchecker::SparseDtmcRegionModelChecker::CoefficientType getConstantPart(storm::RationalFunction const& function){ + return function.constantPart(); + } + template<> bool functionIsLinear(storm::RationalFunction const& function){ // Note: At this moment there is no function in carl for rationalFunctions. diff --git a/src/utility/regions.h b/src/utility/regions.h index ef7a62469..e9e338c3d 100644 --- a/src/utility/regions.h +++ b/src/utility/regions.h @@ -114,9 +114,19 @@ namespace storm { template std::string getVariableName(VariableType variable); + /* + * evaluates the given function at the given point and returns the result + */ template typename storm::modelchecker::SparseDtmcRegionModelChecker::CoefficientType evaluateFunction(ParametricType const& function, std::map::VariableType, typename storm::modelchecker::SparseDtmcRegionModelChecker::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 storm::modelchecker::SparseDtmcRegionModelChecker::CoefficientType getConstantPart(ParametricType const& function); + /*! * Returns true if the function is rational. Note that the function might be simplified. */