diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index 08f09620f..6958f6ac2 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -15,38 +15,53 @@ namespace storm { template - SparseDtmcRegionModelChecker::SamplingModel::SamplingModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : mapping(), evaluationTable(), computeRewards(computeRewards){ + SparseDtmcRegionModelChecker::SamplingModel::SamplingModel(storm::models::sparse::Dtmc const& parametricModel, bool computeRewards) : probabilityMapping(), stateRewardMapping(), evaluationTable(), computeRewards(computeRewards){ // 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; + std::size_t numOfNonConstProbEntries=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; + ++numOfNonConstProbEntries; } matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry); dummyEntry=storm::utility::zero(); } } + this->probabilityMapping.reserve(numOfNonConstProbEntries); - //Obtain other model ingredients and the approximation model itself + //Obtain other model ingredients and the sampling model itself storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling()); - boost::optional> noStateRewards; + boost::optional> stateRewards; + if(this->computeRewards){ + std::size_t numOfNonConstRewEntries=0; + std::vector stateRewardsAsVector(parametricModel.getNumberOfStates()); + for(std::size_t state=0; stateparametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ + stateRewardsAsVector[state] = storm::utility::regions::convertNumber(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state])); + } else { + stateRewardsAsVector[state] = storm::utility::zero(); + functionSet.insert(parametricModel.getStateRewardVector()[state]); + ++numOfNonConstRewEntries; + } + } + stateRewards=std::move(stateRewardsAsVector); + this->stateRewardMapping.reserve(numOfNonConstRewEntries); + } boost::optional> noTransitionRewards; boost::optional>> noChoiceLabeling; - this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(noStateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + this->model=std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); - //Get the evaluation table. Note that it remains sorted due to the fact that functionSet is sorted + //Get the evaluation table. Note that it remains sorted due to the fact that probFunctionSet 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); + //Fill in the entries for the mappings auto samEntry = this->model->getTransitionMatrix().begin(); for(auto const& parEntry : parametricModel.getTransitionMatrix()){ if(this->parametricTypeComparator.isConstant(parEntry.getValue())){ @@ -56,10 +71,20 @@ namespace storm { std::pair searchedPair(parEntry.getValue(), storm::utility::zero()); auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair); STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly"); - mapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); + probabilityMapping.emplace_back(std::make_pair(&(tableIt->second), samEntry)); } ++samEntry; } + if(this->computeRewards){ + for(std::size_t state=0; stateparametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){ + std::pair searchedPair(parametricModel.getStateRewardVector()[state], storm::utility::zero()); + auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair); + STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly"); + stateRewardMapping.emplace_back(std::make_pair(&(tableIt->second), &(this->model->getStateRewardVector()[state]))); + } + } + } } template @@ -79,10 +104,15 @@ namespace storm { 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){ + //write the instantiated values to the matrix according to the mappings + for(auto& mappingPair : this->probabilityMapping){ mappingPair.second->setValue(*(mappingPair.first)); } + if(this->computeRewards){ + for(auto& mappingPair : this->stateRewardMapping){ + *mappingPair.second=*mappingPair.first; + } + } } template diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index 768e82d70..f2ffb0c7d 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -51,7 +51,8 @@ namespace storm { //Vector has one entry for every (non-constant) matrix entry. //pair.first points to an entry in the evaluation table, //pair.second is an iterator to the corresponding matrix entry - std::vector::iterator>> mapping; + std::vector::iterator>> probabilityMapping; + std::vector> stateRewardMapping; //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. diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp index b6305a43a..d7af4d61e 100644 --- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp +++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp @@ -71,49 +71,27 @@ namespace storm { template void SparseDtmcRegionModelChecker::specifyFormula(std::shared_ptr formula) { - std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now(); + std::chrono::high_resolution_clock::time_point timeSpecifyFormulaStart = std::chrono::high_resolution_clock::now(); STORM_LOG_THROW(this->canHandle(*formula), storm::exceptions::InvalidArgumentException, "Tried to specify a formula that can not be handled."); - this->hasOnlyLinearFunctions=false; this->isResultConstant=false; + this->isApproximationApplicable=false; this->smtSolver=nullptr; this->approximationModel=nullptr; this->samplingModel=nullptr; this->reachabilityFunction=nullptr; - //Get bounds, comparison type, target states, .. //Note: canHandle already ensures that the formula has the right shape. this->specifiedFormula = formula; - std::unique_ptr targetStatesResultPtr; - if (this->specifiedFormula->isProbabilityOperatorFormula()) { - this->computeRewards=false; - storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = this->specifiedFormula->asProbabilityOperatorFormula(); - this->specifiedFormulaCompType=probabilityOperatorFormula.getComparisonType(); - this->specifiedFormulaBound=probabilityOperatorFormula.getBound(); - targetStatesResultPtr=this->eliminationModelChecker.check(probabilityOperatorFormula.getSubformula().asEventuallyFormula().getSubformula()); - } - else if (this->specifiedFormula->isRewardOperatorFormula()) { - this->computeRewards=true; - storm::logic::RewardOperatorFormula const& rewardOperatorFormula = this->specifiedFormula->asRewardOperatorFormula(); - this->specifiedFormulaCompType=rewardOperatorFormula.getComparisonType(); - this->specifiedFormulaBound=rewardOperatorFormula.getBound(); - targetStatesResultPtr=this->eliminationModelChecker.check(rewardOperatorFormula.getSubformula().asReachabilityRewardFormula().getSubformula()); - } - else { - STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The specified property " << formula << "is not supported"); - } - storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector(); - - // get a more simple model with a single target state, a single sink state and where states with constant outgoing transitions have been removed - // Note: also checks for linear functions and a constant result - computeSimplifiedModel(targetStates); + // set some information regarding the formula and model. Also computes a more simple version of the model + preprocess(); if(!this->isResultConstant){ //now create the model used for Approximation if(storm::settings::regionSettings().doApprox()){ initializeApproximationModel(*this->simplifiedModel); } //now create the model used for Sampling - if(storm::settings::regionSettings().doSample() || (storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){ + if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE || (storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){ initializeSamplingModel(*this->simplifiedModel); } //Check if the reachability function needs to be computed @@ -123,8 +101,8 @@ namespace storm { } } //some information for statistics... - std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now(); - this->timePreprocessing= timePreprocessingEnd - timePreprocessingStart; + std::chrono::high_resolution_clock::time_point timeSpecifyFormulaEnd = std::chrono::high_resolution_clock::now(); + this->timeSpecifyFormula= timeSpecifyFormulaEnd - timeSpecifyFormulaStart; this->numOfCheckedRegions=0; this->numOfRegionsSolvedThroughSampling=0; this->numOfRegionsSolvedThroughApproximation=0; @@ -140,84 +118,127 @@ namespace storm { } template - void SparseDtmcRegionModelChecker::computeSimplifiedModel(storm::storage::BitVector const& targetStates){ - std::chrono::high_resolution_clock::time_point timeComputeSimplifiedModelStart = std::chrono::high_resolution_clock::now(); - //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; - //if the target states are not reached with probability 1, then the expected reward is defined as infinity - if (this->computeRewards && !this->model.getInitialStates().isSubsetOf(statesWithProbability1)){ - STORM_LOG_WARN("The expected reward of the initial state is constant (infinity)"); - this->reachabilityFunction = std::make_shared(storm::utility::infinity()); - this->isResultConstant=true; - this->hasOnlyLinearFunctions=true; - return; //nothing else to do... + void SparseDtmcRegionModelChecker::preprocess(){ + std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now(); + STORM_LOG_THROW(this->model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidArgumentException, "Input model is required to have exactly one initial state."); + + //first some preprocessing depending on the type of the considered formula + storm::storage::BitVector maybeStates, targetStates; + boost::optional> stateRewards; + if (this->specifiedFormula->isProbabilityOperatorFormula()) { + preprocessForProbabilities(maybeStates, targetStates); } - 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->computeRewards && this->model.getInitialStates().isDisjointFrom(maybeStates)) { - STORM_LOG_WARN("The probability of the initial state is constant (0 or 1)"); - this->reachabilityFunction = std::make_shared(statesWithProbability0.get(*(this->model.getInitialStates().begin())) ? storm::utility::zero() : storm::utility::one()); - this->isResultConstant=true; - this->hasOnlyLinearFunctions=true; - return; //nothing else to do... + else if (this->specifiedFormula->isRewardOperatorFormula()) { + std::vector stateRewardsAsVector; + preprocessForRewards(maybeStates, targetStates, stateRewardsAsVector); + stateRewards=std::move(stateRewardsAsVector); + } + else { + STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The specified property " << this->specifiedFormula << "is not supported"); } // 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); + storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->model.getTransitionMatrix(), this->model.getInitialStates(), maybeStates, targetStates); // 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); + // Create a vector for the probabilities to go to a target state in one step. + std::vector oneStepProbabilities = this->model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, targetStates); // 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. + // eliminate all states with only constant outgoing transitions (and possibly rewards) // Convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. auto flexibleTransitions = this->eliminationModelChecker.getFlexibleSparseMatrix(submatrix); auto 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. - subsystem.set(initialState, false); - this->hasOnlyLinearFunctions=true; - bool allReachableFunctionsAreConstant=true; - boost::optional> missingStateRewards; - for (auto const& state : subsystem) { - bool stateHasOnlyConstantOutgoingTransitions=true; - for(auto const& entry : submatrix.getRow(state)){ + storm::storage::BitVector subsystem(submatrix.getRowCount(), true); + //The states that we consider to eliminate + storm::storage::BitVector considerToEliminate(submatrix.getRowCount(), true); + considerToEliminate.set(initialState, false); + this->isApproximationApplicable=true; + this->isResultConstant=true; + for (auto const& state : considerToEliminate) { + bool eliminateThisState=true; + for(auto const& entry : flexibleTransitions.getRow(state)){ if(!this->parametricTypeComparator.isConstant(entry.getValue())){ - allReachableFunctionsAreConstant=false; - stateHasOnlyConstantOutgoingTransitions=false; - this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear(entry.getValue()); + this->isResultConstant=false; + eliminateThisState=false; + this->isApproximationApplicable &= storm::utility::regions::functionIsLinear(entry.getValue()); } } if(!this->parametricTypeComparator.isConstant(oneStepProbabilities[state])){ - allReachableFunctionsAreConstant=false; - stateHasOnlyConstantOutgoingTransitions=false; - this->hasOnlyLinearFunctions &= storm::utility::regions::functionIsLinear(oneStepProbabilities[state]); + this->isResultConstant=false; + eliminateThisState=false; + this->isApproximationApplicable &= storm::utility::regions::functionIsLinear(oneStepProbabilities[state]); + } + if(this->computeRewards && eliminateThisState && !this->parametricTypeComparator.isConstant(stateRewards.get()[state])){ + //Note: The state reward does not need to be constant but we need to make sure that + //no parameter of this reward function occurs as a parameter in the probability functions of the predecessors + // (otherwise, more complex functions might occur in our simplified model) + // TODO: test if we should really remove these states (the resulting reward functions are less simple this way) + std::set probVars; + for(auto const& predecessor : flexibleBackwardTransitions.getRow(state)){ + for(auto const& predecessorTransition : flexibleTransitions.getRow(predecessor.getColumn())){ + storm::utility::regions::gatherOccurringVariables(predecessorTransition.getValue(), probVars); + } + } + std::set rewardVars; + storm::utility::regions::gatherOccurringVariables(stateRewards.get()[state], rewardVars); + for(auto const& rewardVar : rewardVars){ + if(probVars.find(rewardVar)!=probVars.end()){ + eliminateThisState=false; + break; + } + } } - if(stateHasOnlyConstantOutgoingTransitions){ - this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); + if(eliminateThisState){ + this->eliminationModelChecker.eliminateState(flexibleTransitions, oneStepProbabilities, state, flexibleBackwardTransitions, stateRewards); subsystem.set(state,false); } } - subsystem.set(initialState, true); - STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl); - - if(allReachableFunctionsAreConstant){ + STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions."); + if(this->isResultConstant){ //Check if this is also the case for the initial state - for(auto const& entry : submatrix.getRow(initialState)){ - allReachableFunctionsAreConstant&=this->parametricTypeComparator.isConstant(entry.getValue()); + for(auto const& entry : flexibleTransitions.getRow(initialState)){ + this->isResultConstant&=this->parametricTypeComparator.isConstant(entry.getValue()); + } + this->isResultConstant&=this->parametricTypeComparator.isConstant(oneStepProbabilities[initialState]); + } + if(this->computeRewards && (this->isApproximationApplicable || this->isResultConstant)){ + //We will need to check whether this is also the case for the reward functions. + for(auto const& state : maybeStates){ + std::set rewardVars; + if(this->model.hasStateRewards() && !this->parametricTypeComparator.isConstant(this->model.getStateRewardVector()[state])){ + this->isResultConstant=false; + this->isApproximationApplicable &= storm::utility::regions::functionIsLinear(this->model.getStateRewardVector()[state]); + storm::utility::regions::gatherOccurringVariables(stateRewards.get()[state], rewardVars); + } + if(this->model.hasTransitionRewards()){ + for(auto const& entry : this->model.getTransitionRewardMatrix().getRow(state)) { + if(!this->parametricTypeComparator.isConstant(entry.getValue())){ + this->isResultConstant=false; + this->isApproximationApplicable &= storm::utility::regions::functionIsLinear(entry.getValue()); + storm::utility::regions::gatherOccurringVariables(entry.getValue(), rewardVars); + } + } + } + if(!rewardVars.empty()){ + std::set probVars; + for(auto const& entry: this->model.getTransitionMatrix().getRow(state)){ + storm::utility::regions::gatherOccurringVariables(entry.getValue(), probVars); + } + for(auto const& rewardVar : rewardVars){ + if(probVars.find(rewardVar)!=probVars.end()){ + this->isApproximationApplicable=false; + break; + } + } + } } - allReachableFunctionsAreConstant&=this->parametricTypeComparator.isConstant(oneStepProbabilities[initialState]); - // Set the flag accordingly - this->isResultConstant=allReachableFunctionsAreConstant; - STORM_LOG_WARN_COND(!this->isResultConstant, "For the given property, the reachability probability is constant, i.e., independent of the region"); } + STORM_LOG_WARN_COND(!this->isResultConstant, "For the given property, the reachability Value is constant, i.e., independent of the region"); + //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. @@ -234,21 +255,28 @@ namespace storm { //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()); + ParametricType missingProbability=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(); + missingProbability-=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); + if(this->computeRewards){ + // the missing probability always leads to target + if(!this->parametricTypeComparator.isZero(missingProbability)){ + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState,missingProbability); + } + } else{ + //transition to target state + if(!this->parametricTypeComparator.isZero(oneStepProbabilities[oldStateIndex])){ + missingProbability-=oneStepProbabilities[oldStateIndex]; + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, oneStepProbabilities[oldStateIndex]); + } + //transition to sink state + if(!this->parametricTypeComparator.isZero(missingProbability)){ + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, missingProbability); + } } } //add self loops on the additional states (i.e., target and sink) @@ -266,15 +294,95 @@ namespace storm { sinkLabel.set(sinkState, true); labeling.addLabel("sink", std::move(sinkLabel)); // other ingredients - boost::optional> noStateRewards; + if(this->computeRewards){ + storm::utility::vector::selectVectorValues(stateRewards.get(), subsystem, stateRewards.get()); + stateRewards->push_back(storm::utility::zero()); //target state + stateRewards->push_back(storm::utility::zero()); //sink state + } 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)); - std::chrono::high_resolution_clock::time_point timeComputeSimplifiedModelEnd = std::chrono::high_resolution_clock::now(); - this->timeComputeSimplifiedModel = timeComputeSimplifiedModelEnd - timeComputeSimplifiedModelStart; + this->simplifiedModel = std::make_shared>(matrixBuilder.build(), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling)); + std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now(); + this->timePreprocessing = timePreprocessingEnd - timePreprocessingStart; + } + + template + void SparseDtmcRegionModelChecker::preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates) { + this->computeRewards=false; + //Get bounds, comparison type, target states + storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = this->specifiedFormula->asProbabilityOperatorFormula(); + this->specifiedFormulaCompType=probabilityOperatorFormula.getComparisonType(); + this->specifiedFormulaBound=probabilityOperatorFormula.getBound(); + std::unique_ptr targetStatesResultPtr=this->eliminationModelChecker.check(probabilityOperatorFormula.getSubformula().asEventuallyFormula().getSubformula()); + targetStates = std::move(targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector()); + + //maybeStates: 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); + maybeStates = ~(statesWithProbability01.first | statesWithProbability01.second); + // If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction. + storm::storage::sparse::state_type initialState = *this->model.getInitialStates().begin(); + if (!maybeStates.get(initialState)) { + STORM_LOG_WARN("The probability of the initial state is constant (zero or one)"); + this->reachabilityFunction = std::make_shared(statesWithProbability01.first.get(initialState) ? storm::utility::zero() : storm::utility::one()); + this->isResultConstant=true; + return; //nothing else to do... + } + //extend target states + targetStates=statesWithProbability01.second; + } + + + template + void SparseDtmcRegionModelChecker::preprocessForRewards(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, std::vector& stateRewards) { + this->computeRewards=true; + STORM_LOG_THROW(this->model.hasStateRewards() || this->model.hasTransitionRewards(), storm::exceptions::InvalidArgumentException, "Input model does not have a reward model."); + //Get bounds, comparison type, target states + storm::logic::RewardOperatorFormula const& rewardOperatorFormula = this->specifiedFormula->asRewardOperatorFormula(); + this->specifiedFormulaCompType=rewardOperatorFormula.getComparisonType(); + this->specifiedFormulaBound=rewardOperatorFormula.getBound(); + std::unique_ptr targetStatesResultPtr=this->eliminationModelChecker.check(rewardOperatorFormula.getSubformula().asReachabilityRewardFormula().getSubformula()); + targetStates = std::move(targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector()); + + //maybeStates: Compute the subset of states that has a reachability reward less than infinity. + storm::storage::BitVector statesWithProbability1 = storm::utility::graph::performProb1(this->model.getBackwardTransitions(), storm::storage::BitVector(this->model.getNumberOfStates(), true), targetStates); + maybeStates = ~targetStates & statesWithProbability1; + // If the initial state is known to have 0 reward or an infinite reward value, we can directly set the reachRewardFunction. + storm::storage::sparse::state_type initialState = *this->model.getInitialStates().begin(); + if (!maybeStates.get(initialState)) { + STORM_LOG_WARN("The expected reward of the initial state is constant (infinity or zero)"); + this->reachabilityFunction = std::make_shared(statesWithProbability1.get(initialState) ? storm::utility::zero() : storm::utility::infinity()); + this->isResultConstant=true; + return; //nothing else to do... + } + + stateRewards=std::vector(maybeStates.getNumberOfSetBits()); + if (this->model.hasTransitionRewards()) { + // If a transition-based reward model is available, we initialize the state reward vector to the row sums of the matrix + // that results from taking the pointwise product of the transition probability matrix and the transition reward matrix. + std::vector pointwiseProductRowSumVector = this->model.getTransitionMatrix().getPointwiseProductRowSumVector(this->model.getTransitionRewardMatrix()); + storm::utility::vector::selectVectorValues(stateRewards, maybeStates, pointwiseProductRowSumVector); + + if (this->model.hasStateRewards()) { + // If a state-based reward model is also available, we need to add this vector + // as well. As the state reward vector contains entries not just for the states + // that we still consider (i.e. maybeStates), we need to extract these values + // first. + std::vector subStateRewards(stateRewards.size()); + storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, this->model.getStateRewardVector()); + storm::utility::vector::addVectors(stateRewards, subStateRewards, stateRewards); + } + } else { + // If only a state-based reward model is available, we take this vector as the + // right-hand side. As the state reward vector contains entries not just for the + // states that we still consider (i.e. maybeStates), we need to extract these values + // first. + storm::utility::vector::selectVectorValues(stateRewards, maybeStates, this->model.getStateRewardVector()); + } } + + template void SparseDtmcRegionModelChecker::initializeApproximationModel(storm::models::sparse::Dtmc const& simpleModel) { std::chrono::high_resolution_clock::time_point timeInitApproxModelStart = std::chrono::high_resolution_clock::now(); @@ -311,16 +419,20 @@ namespace storm { //now compute the functions using methods from elimination model checker storm::storage::BitVector newInitialStates = simpleModel.getInitialStates() % maybeStates; storm::storage::BitVector phiStates(simpleModel.getNumberOfStates(), true); - boost::optional> noStateRewards; - + boost::optional> stateRewards; + if(this->computeRewards){ + std::vector stateRewardsAsVector(maybeStates.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(stateRewardsAsVector, maybeStates, simpleModel.getStateRewardVector()); + stateRewards=std::move(stateRewardsAsVector); + } std::vector statePriorities = this->eliminationModelChecker.getStatePriorities(forwardTransitions,backwardTransitions,newInitialStates,oneStepProbabilities); - this->reachabilityFunction=std::make_shared(this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, simpleModel.getStates("target"), noStateRewards, statePriorities)); - /* std::string funcStr = " (/ " + - this->reachProbFunction->nominator().toString(false, true) + " " + - this->reachProbFunction->denominator().toString(false, true) + + this->reachabilityFunction=std::make_shared(this->eliminationModelChecker.computeReachabilityValue(forwardTransitions, oneStepProbabilities, backwardTransitions, newInitialStates , phiStates, simpleModel.getStates("target"), stateRewards, statePriorities)); + /* std::string funcStr = " (/ " + + this->reachabilityFunction->nominator().toString(false, true) + " " + + this->reachabilityFunction->denominator().toString(false, true) + " )"; std::cout << std::endl <<"the resulting reach prob function is " << std::endl << funcStr << std::endl << std::endl;*/ - STORM_LOG_DEBUG("Computed reachProbFunction"); + STORM_LOG_DEBUG("Computed reachabilityFunction"); std::chrono::high_resolution_clock::time_point timeComputeReachabilityFunctionEnd = std::chrono::high_resolution_clock::now(); this->timeComputeReachabilityFunction = timeComputeReachabilityFunctionEnd-timeComputeReachabilityFunctionStart; } @@ -336,8 +448,8 @@ namespace storm { //switches for the different steps. bool done=false; - STORM_LOG_WARN_COND( (!storm::settings::regionSettings().doApprox() || this->hasOnlyLinearFunctions), "the approximation is only correct if the model has only linear functions. As this is not the case, approximation is deactivated"); - bool doApproximation=storm::settings::regionSettings().doApprox() && this->hasOnlyLinearFunctions; + STORM_LOG_WARN_COND( (!storm::settings::regionSettings().doApprox() || this->isApproximationApplicable), "the approximation is only correct if the model has only linear functions. As this is not the case, approximation is deactivated"); + bool doApproximation=storm::settings::regionSettings().doApprox() && this->isApproximationApplicable; bool doSampling=storm::settings::regionSettings().doSample(); bool doFullSmt=storm::settings::regionSettings().doSmt(); @@ -429,7 +541,7 @@ namespace storm { template bool SparseDtmcRegionModelChecker::checkApproximativeValues(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."); + STORM_LOG_THROW(this->isApproximationApplicable, storm::exceptions::UnexpectedException, "Tried to perform approximative method (only applicable if all functions are linear), but there are nonlinear functions."); std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now(); getApproximationModel()->instantiate(region); std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now(); @@ -772,8 +884,8 @@ namespace storm { return; } + std::chrono::milliseconds timeSpecifyFormulaInMilliseconds = std::chrono::duration_cast(this->timeSpecifyFormula); std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast(this->timePreprocessing); - std::chrono::milliseconds timeComputeSimplifiedModelInMilliseconds = std::chrono::duration_cast(this->timeComputeSimplifiedModel); std::chrono::milliseconds timeInitSamplingModelInMilliseconds = std::chrono::duration_cast(this->timeInitSamplingModel); std::chrono::milliseconds timeInitApproxModelInMilliseconds = std::chrono::duration_cast(this->timeInitApproxModel); std::chrono::milliseconds timeComputeReachabilityFunctionInMilliseconds = std::chrono::duration_cast(this->timeComputeReachabilityFunction); @@ -783,7 +895,7 @@ namespace storm { std::chrono::milliseconds timeMDPBuildInMilliseconds = std::chrono::duration_cast(this->timeMDPBuild); std::chrono::milliseconds timeFullSmtInMilliseconds = std::chrono::duration_cast(this->timeFullSmt); - std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ... + std::chrono::high_resolution_clock::duration timeOverall = timeSpecifyFormula + timeCheckRegion; // + ... std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast(timeOverall); uint_fast64_t numOfSolvedRegions= this->numOfRegionsExistsBoth + this->numOfRegionsAllSat + this->numOfRegionsAllViolated; @@ -791,14 +903,14 @@ namespace storm { outstream << std::endl << "Region Model Checker Statistics:" << std::endl; outstream << "-----------------------------------------------" << std::endl; outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl; + outstream << "Formula: " << *this->specifiedFormula << std::endl; if(this->isResultConstant){ outstream << "The requested value is constant (i.e. independent of any parameters)" << std::endl; } else{ outstream << "Simplified model: " << this->simplifiedModel->getNumberOfStates() << " states, " << this->simplifiedModel->getNumberOfTransitions() << " transitions" << std::endl; } - outstream << "Formula: " << *this->specifiedFormula << std::endl; - outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl; + outstream << "Approximation is " << (this->isApproximationApplicable ? "" : "not ") << "applicable" << std::endl; outstream << "Number of checked regions: " << this->numOfCheckedRegions << std::endl; outstream << " Number of solved regions: " << numOfSolvedRegions << "(" << numOfSolvedRegions*100/this->numOfCheckedRegions << "%)" << std::endl; outstream << " AllSat: " << this->numOfRegionsAllSat << "(" << this->numOfRegionsAllSat*100/this->numOfCheckedRegions << "%)" << std::endl; @@ -812,8 +924,8 @@ namespace storm { outstream << std::endl; outstream << "Running times:" << std::endl; outstream << " " << timeOverallInMilliseconds.count() << "ms overall (excluding model parsing, bisimulation (if applied))" << std::endl; - outstream << " " << timePreprocessingInMilliseconds.count() << "ms Preprocessing including... " << std::endl; - outstream << " " << timeComputeSimplifiedModelInMilliseconds.count() << "ms to obtain a simplified model (state elimination of const transitions)" << std::endl; + outstream << " " << timeSpecifyFormulaInMilliseconds.count() << "ms Initialization for the specified formula, including... " << std::endl; + outstream << " " << timePreprocessingInMilliseconds.count() << "ms for Preprocessing (mainly to obtain a simplified model, i.e., state elimination of const transitions)" << std::endl; outstream << " " << timeInitApproxModelInMilliseconds.count() << "ms to initialize the Approximation Model" << std::endl; outstream << " " << timeInitSamplingModelInMilliseconds.count() << "ms to initialize the Sampling Model" << std::endl; outstream << " " << timeComputeReachabilityFunctionInMilliseconds.count() << "ms to compute the reachability function" << std::endl; diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.h b/src/modelchecker/region/SparseDtmcRegionModelChecker.h index d56336974..8ada7b31d 100644 --- a/src/modelchecker/region/SparseDtmcRegionModelChecker.h +++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.h @@ -88,26 +88,54 @@ namespace storm { /*! - * Computes a model with a single target and a single sink state. + * 1. Analyzes the formula (sets this->specifiedFormulaBound, this->specifiedFormulaCompType) + * + * 2. Checks whether the approximation technique is applicable and whether the model checking result is independent of parameters (i.e., constant) + * The flags of This model checker are set accordingly. + * + * 3. Computes a model with a single target and at most one sink state. * Eliminates all states for which the outgoing transitions are constant. + * If rewards are relevant, transition rewards are transformed to state rewards * - * Also checks whether the non constant functions are linear and whether the model checking result is independent of parameters (i.e., constant) - * The flags of This model checker are set accordingly. + * @note this->specifiedFormula has to be set accordingly, before calling this function */ - void computeSimplifiedModel(storm::storage::BitVector const& targetStates); + void preprocess(); + + + /*! + * Does some sanity checks and preprocessing steps on the currently specified model and + * reachability probability formula, i.e., Computes maybeStates and targetStates and sets some flags + * + * @note The returned set of target states also includes states where an 'actual' target state is reached with probability 1 + * + */ + void preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates); + + + /*! + * Does some sanity checks and preprocessing steps on the currently specified model and + * reachability reward formula, i.e., computes some flags, maybeStates, targetStates and + * a new stateReward vector that considers state+transition rewards of the original model. + * @note stateRewards.size = maybeStates.numberOfSetBits + */ + void preprocessForRewards(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, std::vector& stateRewards); /*! * initializes the Approximation Model + * + * @note computeFlagsAndSimplifiedModel should be called before calling this */ void initializeApproximationModel(storm::models::sparse::Dtmc const& simpleModel); /*! * initializes the Sampling Model + * @note computeFlagsAndSimplifiedModel should be called before calling this */ void initializeSamplingModel(storm::models::sparse::Dtmc const& simpleModel); /*! * Computes the reachability function via state elimination + * @note computeFlagsAndSimplifiedModel should be called before calling this */ void computeReachabilityFunction(storm::models::sparse::Dtmc const& simpleModel); @@ -208,7 +236,7 @@ namespace storm { // The function for the reachability probability (or: reachability reward) in the initial state std::shared_ptr reachabilityFunction; // a flag that is true if there are only linear functions at transitions of the model - bool hasOnlyLinearFunctions; + bool isApproximationApplicable; // a flag that is true iff the resulting reachability function is constant bool isResultConstant; // the smt solver that is used to prove properties with the help of the reachabilityFunction @@ -223,8 +251,8 @@ namespace storm { uint_fast64_t numOfRegionsAllSat; uint_fast64_t numOfRegionsAllViolated; + std::chrono::high_resolution_clock::duration timeSpecifyFormula; std::chrono::high_resolution_clock::duration timePreprocessing; - std::chrono::high_resolution_clock::duration timeComputeSimplifiedModel; std::chrono::high_resolution_clock::duration timeInitApproxModel; std::chrono::high_resolution_clock::duration timeInitSamplingModel; std::chrono::high_resolution_clock::duration timeComputeReachabilityFunction; diff --git a/src/models/sparse/Model.cpp b/src/models/sparse/Model.cpp index 506083cbc..3e8d64b2a 100644 --- a/src/models/sparse/Model.cpp +++ b/src/models/sparse/Model.cpp @@ -103,6 +103,11 @@ namespace storm { return stateRewardVector.get(); } + template + std::vector& Model::getStateRewardVector() { + return stateRewardVector.get(); + } + template boost::optional> const& Model::getOptionalStateRewardVector() const { return stateRewardVector; diff --git a/src/models/sparse/Model.h b/src/models/sparse/Model.h index 83173ab63..606bb104a 100644 --- a/src/models/sparse/Model.h +++ b/src/models/sparse/Model.h @@ -168,6 +168,14 @@ namespace storm { */ std::vector const& getStateRewardVector() const; + /*! + * Retrieves a vector representing the state rewards of the model. Note that calling this method is only + * valid if the model has state rewards. + * + * @return A vector representing the state rewards of the model. + */ + std::vector& getStateRewardVector(); + /*! * Retrieves an optional value that contains the state reward vector if there is one. * diff --git a/src/settings/modules/RegionSettings.cpp b/src/settings/modules/RegionSettings.cpp index 9c9d14b09..633754132 100644 --- a/src/settings/modules/RegionSettings.cpp +++ b/src/settings/modules/RegionSettings.cpp @@ -27,7 +27,7 @@ namespace storm { std::vector approxModes = {"off", "guessallsat", "guessallviolated", "testfirst"}; this->addOption(storm::settings::OptionBuilder(moduleName, approxmodeOptionName, true, "Sets whether approximation should be done and whether lower or upper bounds are computed first.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("mode", "The mode, (off, guessallsat, guessallviolated, testfirst). E.g. guessallsat will first try to prove ALLSAT") - .addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(approxModes)).setDefaultValueString("guessallsat").build()).build()); + .addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(approxModes)).setDefaultValueString("testfirst").build()).build()); std::vector sampleModes = {"off", "instantiate", "evaluate"}; this->addOption(storm::settings::OptionBuilder(moduleName, samplemodeOptionName, true, "Sets whether sampling should be done and whether to instantiate a model or compute+evaluate a function.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("mode", "The mode, (off, instantiate, evaluate)")