diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index 6571e45e9..31ee02620 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -16,9 +16,6 @@ #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/InvalidStateException.h" -#include "exceptions/UnexpectedException.h" -#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h" -#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h" namespace storm { namespace modelchecker { @@ -987,111 +984,7 @@ namespace storm { } std::cout << std::endl; } - } - -#ifdef STORM_HAVE_CARL - template<> - std::pair,std::vector>> SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const { - - //Check if the arguments are as expected - STORM_LOG_THROW((std::is_same::value), storm::exceptions::IllegalArgumentException, "Unexpected Type of Coefficients"); - STORM_LOG_THROW(filter.size()==this->getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the filter"); - STORM_LOG_THROW(oneStepProbabilities.empty() || oneStepProbabilities.size()==this->getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the oneStepProbabilities"); - - //get a mapping from old state indices to the new ones - std::vector newStateIndexMap(this->getNumberOfRows(), this->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 : filter){ - newStateIndexMap[state]=newStateIndex; - ++newStateIndex; - } - index_type numStates=filter.getNumberOfSetBits(); - STORM_LOG_ASSERT(newStateIndex==numStates, "unexpected number of new states"); - storm::storage::sparse::state_type targetState =0; - storm::storage::sparse::state_type sinkState=0; - if(!oneStepProbabilities.empty()){ - targetState=numStates; - ++numStates; - } - if(addSinkState){ - sinkState=numStates; - ++numStates; - } - //todo rows (i.e. the first parameter) should be numStates*substitutions.size ? - storm::storage::SparseMatrixBuilder matrixBuilder(0, numStates, 0, true, true, numStates); - std::vector> choiceLabeling; - //fill in the data row by row - index_type currentRow=0; - for(auto const& oldStateIndex : filter){ - matrixBuilder.newRowGroup(currentRow); - for(index_type substitutionIndex=0; substitutionIndex< substitutions.size(); ++substitutionIndex){ - double missingProbability = 1.0; - if(this->getRow(oldStateIndex).empty()){ //just add the selfloop if there is no transition - if(addSelfLoops){ - matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero()); - } - } - else{ - const_iterator entry = this->getRow(oldStateIndex).begin(); - for(; entrygetRow(oldStateIndex).end() && entry->getColumn()getValue().evaluate(substitutions[substitutionIndex])); - missingProbability-=value; - storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; - STORM_LOG_THROW(columngetColumn()!=oldStateIndex){ //maybe add a zero valued diagonal entry - matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero()); - } - for(; entry < this->getRow(oldStateIndex).end(); ++entry){ //insert the rest - double value = cln::double_approx(entry->getValue().evaluate(substitutions[substitutionIndex])); - missingProbability-=value; - storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; - STORM_LOG_THROW(column doubleComperator; - if(addSinkState && !doubleComperator.isZero(missingProbability)){ //transition to sink state - STORM_LOG_ASSERT(missingProbability> -storm::settings::generalSettings().getPrecision(), "The missing probability is negative."); - matrixBuilder.addNextValue(currentRow, sinkState, missingProbability); - } - boost::container::flat_set label; - label.insert(substitutionIndex); - choiceLabeling.emplace_back(label); - ++currentRow; - } - } - //finally, add self loops on the additional states (i.e., target and sink) - boost::container::flat_set labelAll; - labelAll.insert(substitutions.size()); - if (!oneStepProbabilities.empty()){ - matrixBuilder.newRowGroup(currentRow); - matrixBuilder.addNextValue(currentRow, targetState, storm::utility::one()); - choiceLabeling.emplace_back(labelAll); - ++currentRow; - } - - if (addSinkState){ - matrixBuilder.newRowGroup(currentRow); - matrixBuilder.addNextValue(currentRow, sinkState, storm::utility::one()); - choiceLabeling.emplace_back(labelAll); - ++currentRow; - } - - return std::pair, std::vector>>(matrixBuilder.build(), std::move(choiceLabeling)); - } - - template - std::pair,std::vector>> SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const{ STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Instantiation of flexible matrix is not supported for this type"); - } -#endif - + } template typename SparseDtmcEliminationModelChecker::FlexibleSparseMatrix SparseDtmcEliminationModelChecker::getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne) { @@ -1120,388 +1013,11 @@ namespace storm { return flexibleMatrix; } -#ifdef STORM_HAVE_CARL - - template<> - void SparseDtmcEliminationModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleSparseMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleSparseMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities){ - - if(true){ // eliminate all states with constant outgoing transitions - storm::storage::BitVector statesToEliminate = ~initialstates; - //todo: ordering of states important? - boost::optional> missingStateRewards; - for (auto const& state : statesToEliminate) { - bool onlyConstantOutgoingTransitions=true; - for(auto const& entry : flexibleMatrix.getRow(state)){ - if(!entry.getValue().isConstant()){ - onlyConstantOutgoingTransitions=false; - break; - } - } - if(onlyConstantOutgoingTransitions){ - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); - subsystem.set(state,false); - } - } - - //Note: we could also "eliminate" the initial state to get rid of its selfloop - } - else if(false){ //eliminate all states with standard state elimination - boost::optional> missingStateRewards; - - storm::storage::BitVector statesToEliminate = ~initialstates; - std::vector states(statesToEliminate.begin(), statesToEliminate.end()); - - if (statePriorities) { - std::sort(states.begin(), states.end(), [&statePriorities] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statePriorities.get()[a] < statePriorities.get()[b]; }); - } - - STORM_LOG_DEBUG("Eliminating " << states.size() << " states using the state elimination technique." << std::endl); - for (auto const& state : states) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); - } - subsystem=~statesToEliminate; - - } - else if(false){ //hybrid method - boost::optional> missingStateRewards; - storm::storage::BitVector statesToEliminate = ~initialstates; - uint_fast64_t maximalDepth =0; - std::vector entryStateQueue; - STORM_LOG_DEBUG("Eliminating " << statesToEliminate.size() << " states using the hybrid elimination technique." << std::endl); - maximalDepth = treatScc(flexibleMatrix, oneStepProbabilities, initialstates, statesToEliminate, forwardTransitions, flexibleBackwardTransitions, 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) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); - } - } - subsystem=~statesToEliminate; - } - std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << "states." << std::endl; - STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " states." << std::endl); - } - - template - void SparseDtmcEliminationModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleSparseMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleSparseMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities){ - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "elimination of states not suported for this type"); - } - - - template<> - void SparseDtmcEliminationModelChecker::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities){ - carl::VariablePool& varPool = carl::VariablePool::getInstance(); - - //first add a state variable for every state in the subsystem, providing that such a variable does not already exist. - for (storm::storage::sparse::state_type state : subsystem){ - if(stateProbVars[state].isZero()){ //variable does not exist yet - storm::Variable stateVar = varPool.getFreshVariable("p_" + std::to_string(state)); - std::shared_ptr>> cache(new carl::Cache>()); - storm::RationalFunction::PolyType stateVarAsPoly(storm::RationalFunction::PolyType::PolyType(stateVar), cache); - - //each variable is in the interval [0,1] - solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::GEQ, storm::RationalFunction(0)); - solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::LEQ, storm::RationalFunction(1)); - stateProbVars[state] = stateVarAsPoly; - } - } - - //now lets add the actual transitions - for (storm::storage::sparse::state_type state : subsystem){ - storm::RationalFunction reachProbability(oneStepProbabilities[state]); - for(auto const& transition : flexibleMatrix.getRow(state)){ - reachProbability += transition.getValue() * stateProbVars[transition.getColumn()]; - } - //Todo: depending on the objective (i.e. the formlua) it suffices to use LEQ or GEQ here... maybe this is faster? - solver.add(storm::RationalFunction(stateProbVars[state]), storm::CompareRelation::EQ, reachProbability); - } - } - - template - void SparseDtmcEliminationModelChecker::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities){ - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "SMT formulation is not supported for this type"); - } - - template<> - void SparseDtmcEliminationModelChecker::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compType){ - //We are going to build a new (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter - - //todo invent something better to obtain the substitutions. - //this only works as long as there is only one parameter per state, - // also: check whether the terms are linear/monotone(?) - - STORM_LOG_WARN("the probability restriction is experimental.. only works on linear terms with one parameter per state"); - storm::storage::sparse::state_type const numOfStates=subsystem.getNumberOfSetBits() + 2; //subsystem + target state + sink state - storm::models::sparse::StateLabeling stateLabeling(numOfStates); - stateLabeling.addLabel("init", storm::storage::BitVector(numOfStates, true)); - storm::storage::BitVector targetLabel(numOfStates, false); - targetLabel.set(numOfStates-2, true); - stateLabeling.addLabel("target", std::move(targetLabel)); - storm::storage::BitVector sinkLabel(numOfStates, false); - sinkLabel.set(numOfStates-1, true); - stateLabeling.addLabel("sink", std::move(sinkLabel)); - - std::map substitutionLB; - for(auto const& parRegion : regions){ - substitutionLB.insert(std::pair(parRegion.variable, parRegion.lowerBound)); - } - std::map substitutionUB; - for(auto const& parRegion : regions){ - substitutionUB.insert(std::pair(parRegion.variable, parRegion.upperBound)); - } - std::vector> substitutions; - substitutions.push_back(substitutionLB); - substitutions.push_back(substitutionUB); - - std::pair,std::vector>> instantiation = flexibleMatrix.instantiateAsDouble(substitutions, subsystem, true, oneStepProbabilities, true); - boost::optional> noStateRewards; - boost::optional> noTransitionRewards; - storm::models::sparse::Mdp mdp(instantiation.first, std::move(stateLabeling),noStateRewards,noTransitionRewards,instantiation.second); - - //we need the correct optimalityType for model checking as well as the correct relation for smt solving - storm::logic::OptimalityType opType; - storm::CompareRelation boundRelation; - switch (compType){ - case storm::logic::ComparisonType::Greater: - opType=storm::logic::OptimalityType::Minimize; - boundRelation=storm::CompareRelation::GEQ; - break; - case storm::logic::ComparisonType::GreaterEqual: - opType=storm::logic::OptimalityType::Minimize; - boundRelation=storm::CompareRelation::GEQ; - break; - case storm::logic::ComparisonType::Less: - opType=storm::logic::OptimalityType::Maximize; - boundRelation=storm::CompareRelation::LEQ; - break; - case storm::logic::ComparisonType::LessEqual: - opType=storm::logic::OptimalityType::Maximize; - boundRelation=storm::CompareRelation::LEQ; - break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); - } - - //perform model checking on the mdp - std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); - storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); - storm::modelchecker::SparseMdpPrctlModelChecker modelChecker(mdp); - std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType); - std::vector resultVector = resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); - - //todo this is experimental.. - if (true){ - std::cout << "the matrix has " << mdp.getTransitionMatrix().getRowGroupCount() << "row groups, " << mdp.getTransitionMatrix().getRowCount() << " rows, and " << mdp.getTransitionMatrix().getColumnCount() << "columns" << std::endl; - boost::container::flat_set< uint_fast64_t> lbChoiceLabel; - lbChoiceLabel.insert(1); - lbChoiceLabel.insert(2); - storm::models::sparse::Dtmc dtmc(mdp.restrictChoiceLabels(lbChoiceLabel).getTransitionMatrix(), std::move(stateLabeling)); - //modelchecking on dtmc - storm::modelchecker::SparseDtmcPrctlModelChecker dtmcModelChecker(dtmc); - std::unique_ptr resultPtrDtmc = dtmcModelChecker.computeEventuallyProbabilities(eventuallyFormula); - std::vector resultVectorDtmc = resultPtrDtmc->asExplicitQuantitativeCheckResult().getValueVector(); - std::cout << "dtmc result with lower bounds:" << resultVectorDtmc[0]; - std::cout << "mdp result:" << resultVector[0]; - } - - - //formulate constraints for the solver - uint_fast64_t boundDenominator = 1.0/storm::settings::generalSettings().getPrecision(); //we need to approx. the obtained bounds as rational numbers - storm::storage::sparse::state_type subsystemState=0; //the subsystem uses other state indices - for(storm::storage::sparse::state_type state : subsystem){ - uint_fast64_t boundNumerator = resultVector[subsystemState]*boundDenominator; - storm::RationalFunction bound(boundNumerator); - bound = bound/boundDenominator; - //Todo: non-exact values might be problematic here... - solver.add(storm::RationalFunction(stateProbVars[state]), boundRelation, bound); - ++subsystemState; - } - } - - template - void SparseDtmcEliminationModelChecker::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compType){ - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "restricting Probability Variables is not supported for this type"); - } - - template<> - bool SparseDtmcEliminationModelChecker::checkRegion(storm::logic::Formula const& formula, std::vector::ParameterRegion> parameterRegions){ - //Note: this is an 'experimental' implementation - - std::chrono::high_resolution_clock::time_point timeStart = std::chrono::high_resolution_clock::now(); - - //Start with some preprocessing (inspired by computeUntilProbabilities...) - //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ]) - //get the (sub)formulae and the vector of target states - STORM_LOG_THROW(formula.isStateFormula(), storm::exceptions::IllegalArgumentException, "expected a stateFormula"); - STORM_LOG_THROW(formula.asStateFormula().isProbabilityOperatorFormula(), storm::exceptions::IllegalArgumentException, "expected a probabilityOperatorFormula"); - storm::logic::ProbabilityOperatorFormula const& probOpForm=formula.asStateFormula().asProbabilityOperatorFormula(); - STORM_LOG_THROW(probOpForm.hasBound(), storm::exceptions::IllegalArgumentException, "The formula has no bound"); - STORM_LOG_THROW(probOpForm.getSubformula().asPathFormula().isEventuallyFormula(), storm::exceptions::IllegalArgumentException, "expected an eventually subformula"); - storm::logic::EventuallyFormula const& eventuallyFormula = probOpForm.getSubformula().asPathFormula().asEventuallyFormula(); - std::unique_ptr targetStatesResultPtr = this->check(eventuallyFormula.getSubformula()); - storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector(); - // Do some sanity checks to establish some required properties. - STORM_LOG_THROW(model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); - storm::storage::sparse::state_type initialState = *model.getInitialStates().begin(); - // 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 return the result. - if (model.getInitialStates().isDisjointFrom(maybeStates)) { - STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step."); - double res= statesWithProbability0.get(*model.getInitialStates().begin()) ? 0.0 : 1.0; - switch (probOpForm.getComparisonType()){ - case storm::logic::ComparisonType::Greater: - return (res > probOpForm.getBound()); - case storm::logic::ComparisonType::GreaterEqual: - return (res >= probOpForm.getBound()); - case storm::logic::ComparisonType::Less: - return (res < probOpForm.getBound()); - case storm::logic::ComparisonType::LessEqual: - return (res <= probOpForm.getBound()); - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula 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(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. - std::vector oneStepProbabilities = model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); - // Determine the set of initial states of the sub-model. - storm::storage::BitVector newInitialStates = model.getInitialStates() % maybeStates; - // 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(); - - std::vector statePriorities = getStatePriorities(submatrix, submatrixTransposed, newInitialStates, oneStepProbabilities); - - // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. - FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(submatrix); - FlexibleSparseMatrix flexibleBackwardTransitions = getFlexibleSparseMatrix(submatrixTransposed, true); - - std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now(); - - // 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); - eliminateStates(subsystem, flexibleMatrix, oneStepProbabilities, flexibleBackwardTransitions, newInitialStates, submatrix, statePriorities); - - std::chrono::high_resolution_clock::time_point timeStateElemEnd = std::chrono::high_resolution_clock::now(); - - // SMT formulation of resulting pdtmc - storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions - storm::solver::Smt2SmtSolver solver(manager, true); - // we will introduce a variable for every state which encodes the probability to reach a target state from this state. - // we will store them as polynomials to easily use operations with rational functions - std::vector stateProbVars(subsystem.size(), storm::RationalFunction::PolyType(0)); - // todo maybe introduce the parameters already at this point? - formulateModelWithSMT(solver, stateProbVars, subsystem, flexibleMatrix, oneStepProbabilities); - - //the property should be satisfied in the initial state for all parameters. - //this is equivalent to: - //the negation of the property should not be satisfied for some parameter valuation. - //Hence, we flip the comparison relation and later check whether all the constraints are unsat. - storm::CompareRelation propertyCompRel; - switch (probOpForm.getComparisonType()){ - case storm::logic::ComparisonType::Greater: - propertyCompRel=storm::CompareRelation::LEQ; - break; - case storm::logic::ComparisonType::GreaterEqual: - propertyCompRel=storm::CompareRelation::LT; - break; - case storm::logic::ComparisonType::Less: - propertyCompRel=storm::CompareRelation::GEQ; - break; - case storm::logic::ComparisonType::LessEqual: - propertyCompRel=storm::CompareRelation::GT; - break; - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); - } - uint_fast64_t thresholdDenominator = 1.0/storm::settings::generalSettings().getPrecision(); - uint_fast64_t thresholdNumerator = probOpForm.getBound()*thresholdDenominator; - storm::RationalFunction threshold(thresholdNumerator); - threshold = threshold / thresholdDenominator; - solver.add(storm::RationalFunction(stateProbVars[*newInitialStates.begin()]), propertyCompRel, threshold); - - //the bounds for the parameters - solver.push(); - for(auto param : parameterRegions){ - storm::RawPolynomial lB(param.variable); - lB -= param.lowerBound; - solver.add(carl::Constraint(lB,storm::CompareRelation::GEQ)); - storm::RawPolynomial uB(param.variable); - uB -= param.upperBound; - solver.add(carl::Constraint(uB,storm::CompareRelation::LEQ)); - } - - std::chrono::high_resolution_clock::time_point timeSmtFormulationEnd = std::chrono::high_resolution_clock::now(); - - // find further restriction on probabilities - restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions, storm::logic::ComparisonType::Less); //probOpForm.getComparisonType()); - restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions, storm::logic::ComparisonType::Greater); - - std::chrono::high_resolution_clock::time_point timeRestrictingEnd = std::chrono::high_resolution_clock::now(); - - std::cout << "start solving ..." << std::endl; - bool result; - switch (solver.check()){ - case storm::solver::SmtSolver::CheckResult::Sat: - std::cout << "sat!" << std::endl; - result=false; - break; - case storm::solver::SmtSolver::CheckResult::Unsat: - std::cout << "unsat!" << std::endl; - result=true; - break; - case storm::solver::SmtSolver::CheckResult::Unknown: - std::cout << "unknown!" << std::endl; - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (Check-result: Unknown)") - result=false; - break; - default: - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (unexpected check-result)") - result=false; - } - - std::chrono::high_resolution_clock::time_point timeSolvingEnd = std::chrono::high_resolution_clock::now(); - - std::chrono::high_resolution_clock::duration timePreprocessing = timePreprocessingEnd - timeStart; - std::chrono::high_resolution_clock::duration timeStateElem = timeStateElemEnd - timePreprocessingEnd; - std::chrono::high_resolution_clock::duration timeSmtFormulation = timeSmtFormulationEnd - timeStateElemEnd; - std::chrono::high_resolution_clock::duration timeRestricting = timeRestrictingEnd - timeSmtFormulationEnd; - std::chrono::high_resolution_clock::duration timeSolving = timeSolvingEnd- timeRestrictingEnd; - std::chrono::high_resolution_clock::duration timeOverall = timeSolvingEnd - timeStart; - std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast(timePreprocessing); - std::chrono::milliseconds timeStateElemInMilliseconds = std::chrono::duration_cast(timeStateElem); - std::chrono::milliseconds timeSmtFormulationInMilliseconds = std::chrono::duration_cast(timeSmtFormulation); - std::chrono::milliseconds timeRestrictingInMilliseconds = std::chrono::duration_cast(timeRestricting); - std::chrono::milliseconds timeSolvingInMilliseconds = std::chrono::duration_cast(timeSolving); - std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast(timeOverall); - STORM_PRINT_AND_LOG(std::endl << "required time: " << timeOverallInMilliseconds.count() << "ms. Time Breakdown:" << std::endl); - STORM_PRINT_AND_LOG(" * " << timePreprocessingInMilliseconds.count() << "ms for Preprocessing" << std::endl); - STORM_PRINT_AND_LOG(" * " << timeStateElemInMilliseconds.count() << "ms for StateElemination" << std::endl); - STORM_PRINT_AND_LOG(" * " << timeSmtFormulationInMilliseconds.count() << "ms for SmtFormulation" << std::endl); - STORM_PRINT_AND_LOG(" * " << timeRestrictingInMilliseconds.count() << "ms for Restricting" << std::endl); - STORM_PRINT_AND_LOG(" * " << timeSolvingInMilliseconds.count() << "ms for Solving" << std::endl); - - return result; - } - - template - bool SparseDtmcEliminationModelChecker::checkRegion(storm::logic::Formula const& formula, std::vector::ParameterRegion> parameterRegions){ - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Region check is not supported for this type"); - } -#endif template class SparseDtmcEliminationModelChecker; #ifdef STORM_HAVE_CARL template class SparseDtmcEliminationModelChecker; #endif } // namespace modelchecker -} // namespace storm +} // namespace storm \ No newline at end of file diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h index 0745b3e1b..3601ef46d 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h @@ -3,17 +3,23 @@ #include "src/storage/sparse/StateType.h" #include "src/models/sparse/Dtmc.h" -#include "src/models/sparse/Mdp.h" #include "src/modelchecker/AbstractModelChecker.h" #include "src/utility/constants.h" -#include "src/solver/SmtSolver.h" -#include "src/solver/Smt2SmtSolver.h" +//forward declaration of friend class +namespace storm { + namespace modelchecker { + template + class SparseDtmcRegionModelChecker; + } +} + namespace storm { namespace modelchecker { template class SparseDtmcEliminationModelChecker : public AbstractModelChecker { + template friend class storm::modelchecker::SparseDtmcRegionModelChecker; public: explicit SparseDtmcEliminationModelChecker(storm::models::sparse::Dtmc const& model); @@ -24,22 +30,6 @@ namespace storm { virtual std::unique_ptr computeConditionalProbabilities(storm::logic::ConditionalPathFormula const& pathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr checkBooleanLiteralFormula(storm::logic::BooleanLiteralFormula const& stateFormula) override; virtual std::unique_ptr checkAtomicLabelFormula(storm::logic::AtomicLabelFormula const& stateFormula) override; - -#ifdef STORM_HAVE_CARL - struct ParameterRegion{ - storm::Variable variable; - storm::RationalFunction::CoeffType lowerBound; - storm::RationalFunction::CoeffType upperBound; - }; - - - /*! - * Checks whether the given formula holds for all possible parameters that satisfy the given parameter regions - * ParameterRegions should contain all parameters (not mentioned parameters are assumed to be arbitrary reals) - */ - bool checkRegion(storm::logic::Formula const& formula, std::vector parameterRegions); -#endif - private: class FlexibleSparseMatrix { @@ -71,29 +61,6 @@ namespace storm { */ bool hasSelfLoop(storm::storage::sparse::state_type state) const; -#ifdef STORM_HAVE_CARL - /*! - * Instantiates the matrix, i.e., evaluate the occurring functions according to the given substitutions of the variables. - * If there are multiple substitutions, the matrix will consist of multiple row groups. It can be seen as the transition matrix of an MDP - * Only the rows selected by the given filter are considered. (filter should have size==this->getNumberOfRows()) - * An exception is thrown if there is a transition from a selected state to an unselected state - * If one step probabilities are given, a new state is added which can be considered as target state. - * The "missing" probability can be redirected to a sink state - * By convention, the target state will have index filter.getNumberOfSetBits() and the sink state will be the state with the highest index (so right after the target state) - * - * - * @param substitutions A list of mappings, each assigning a constant value to every variable - * @param filter selects the rows of this flexibleMatrix, that will be considered - * @param addSinkState adds a state with a self loop to which the "missing" probability will lead - * @param oneStepProbabilities if given, a new state is added to which there are transitions for all non-zero entries in this vector - * @param addSelfLoops if set, zero valued selfloops will be added in every row - * - * @return A matrix with constant (double) entries and a choice labeling - */ - std::pair,std::vector>> instantiateAsDouble(std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState=true, std::vector const& oneStepProbabilities=std::vector(), bool addSelfLoops=true) const; - //todo add const keyword -#endif - private: std::vector data; }; @@ -108,13 +75,6 @@ namespace storm { std::vector getStatePriorities(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& transitionMatrixTransposed, storm::storage::BitVector const& initialStates, std::vector const& oneStepProbabilities); - //eliminates some of the states according to different strategies. - void eliminateStates(storm::storage::BitVector& subsystem, FlexibleSparseMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleSparseMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities = {}); - - void formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities); - - void restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleSparseMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compTypeOfProperty); - // The model this model checker is supposed to analyze. storm::models::sparse::Dtmc const& model; @@ -125,4 +85,4 @@ namespace storm { } // namespace modelchecker } // namespace storm -#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCELIMINATIONMODELCHECKER_H_ */ +#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCELIMINATIONMODELCHECKER_H_ */ \ No newline at end of file diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp new file mode 100644 index 000000000..3ac5a0ce9 --- /dev/null +++ b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp @@ -0,0 +1,557 @@ +#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h" + +//TODO remove useless includes + +//#include +#include + +#include "src/adapters/CarlAdapter.h" + +//#include "src/storage/StronglyConnectedComponentDecomposition.h" + +#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" + +#include "src/utility/graph.h" +#include "src/utility/vector.h" +#include "src/utility/macros.h" + +#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h" +#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h" +//#include "modelchecker/reachability/SparseDtmcEliminationModelChecker.h" + +#include "src/exceptions/InvalidPropertyException.h" +#include "src/exceptions/InvalidStateException.h" +#include "exceptions/UnexpectedException.h" + +namespace storm { + namespace modelchecker { + + template + SparseDtmcRegionModelChecker::SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc const& model) : model(model), eliminationModelChecker(model) { + // Intentionally left empty. + } + + template + //Todo adapt + bool SparseDtmcRegionModelChecker::canHandle(storm::logic::Formula const& formula) const { + if (formula.isProbabilityOperatorFormula()) { + storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = formula.asProbabilityOperatorFormula(); + return this->canHandle(probabilityOperatorFormula.getSubformula()); + } else if (formula.isRewardOperatorFormula()) { + storm::logic::RewardOperatorFormula const& rewardOperatorFormula = formula.asRewardOperatorFormula(); + return this->canHandle(rewardOperatorFormula.getSubformula()); + } else if (formula.isUntilFormula() || formula.isEventuallyFormula()) { + if (formula.isUntilFormula()) { + storm::logic::UntilFormula const& untilFormula = formula.asUntilFormula(); + if (untilFormula.getLeftSubformula().isPropositionalFormula() && untilFormula.getRightSubformula().isPropositionalFormula()) { + return true; + } + } else if (formula.isEventuallyFormula()) { + storm::logic::EventuallyFormula const& eventuallyFormula = formula.asEventuallyFormula(); + if (eventuallyFormula.getSubformula().isPropositionalFormula()) { + return true; + } + } + } else if (formula.isReachabilityRewardFormula()) { + storm::logic::ReachabilityRewardFormula reachabilityRewardFormula = formula.asReachabilityRewardFormula(); + if (reachabilityRewardFormula.getSubformula().isPropositionalFormula()) { + return true; + } + } else if (formula.isConditionalPathFormula()) { + storm::logic::ConditionalPathFormula conditionalPathFormula = formula.asConditionalPathFormula(); + if (conditionalPathFormula.getLeftSubformula().isEventuallyFormula() && conditionalPathFormula.getRightSubformula().isEventuallyFormula()) { + return this->canHandle(conditionalPathFormula.getLeftSubformula()) && this->canHandle(conditionalPathFormula.getRightSubformula()); + } + } else if (formula.isPropositionalFormula()) { + return true; + } + return false; + } + + template<> + std::pair,std::vector>> SparseDtmcRegionModelChecker::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const { + + //Check if the arguments are as expected + STORM_LOG_THROW((std::is_same::value), storm::exceptions::IllegalArgumentException, "Unexpected Type of Coefficients"); + STORM_LOG_THROW(filter.size()==matrix.getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the filter"); + STORM_LOG_THROW(oneStepProbabilities.empty() || oneStepProbabilities.size()==matrix.getNumberOfRows(), storm::exceptions::IllegalArgumentException, "Unexpected size of the oneStepProbabilities"); + + //get a mapping from old state indices to the new ones + std::vector newStateIndexMap(matrix.getNumberOfRows(), matrix.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 : filter){ + newStateIndexMap[state]=newStateIndex; + ++newStateIndex; + } + storm::storage::sparse::state_type numStates=filter.getNumberOfSetBits(); + STORM_LOG_ASSERT(newStateIndex==numStates, "unexpected number of new states"); + storm::storage::sparse::state_type targetState =0; + storm::storage::sparse::state_type sinkState=0; + if(!oneStepProbabilities.empty()){ + targetState=numStates; + ++numStates; + } + if(addSinkState){ + sinkState=numStates; + ++numStates; + } + //todo rows (i.e. the first parameter) should be numStates*substitutions.size ? + 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(auto const& oldStateIndex : filter){ + matrixBuilder.newRowGroup(currentRow); + for(size_t substitutionIndex=0; substitutionIndex< substitutions.size(); ++substitutionIndex){ + double missingProbability = 1.0; + if(matrix.getRow(oldStateIndex).empty()){ //just add the selfloop if there is no transition + if(addSelfLoops){ + matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero()); + } + } + else{ + FlexibleMatrix::const_iterator entry = matrix.getRow(oldStateIndex).begin(); + for(; entrygetColumn()getValue().evaluate(substitutions[substitutionIndex])); + missingProbability-=value; + storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; + STORM_LOG_THROW(columngetColumn()!=oldStateIndex){ //maybe add a zero valued diagonal entry + matrixBuilder.addNextValue(currentRow, newStateIndexMap[oldStateIndex], storm::utility::zero()); + } + for(; entry < matrix.getRow(oldStateIndex).end(); ++entry){ //insert the rest + double value = cln::double_approx(entry->getValue().evaluate(substitutions[substitutionIndex])); + missingProbability-=value; + storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; + STORM_LOG_THROW(column doubleComperator; + if(addSinkState && !doubleComperator.isZero(missingProbability)){ //transition to sink state + STORM_LOG_ASSERT(missingProbability> -storm::settings::generalSettings().getPrecision(), "The missing probability is negative."); + matrixBuilder.addNextValue(currentRow, sinkState, missingProbability); + } + boost::container::flat_set label; + label.insert(substitutionIndex); + choiceLabeling.emplace_back(label); + ++currentRow; + } + } + //finally, add self loops on the additional states (i.e., target and sink) + boost::container::flat_set labelAll; + labelAll.insert(substitutions.size()); + if (!oneStepProbabilities.empty()){ + matrixBuilder.newRowGroup(currentRow); + matrixBuilder.addNextValue(currentRow, targetState, storm::utility::one()); + choiceLabeling.emplace_back(labelAll); + ++currentRow; + } + + if (addSinkState){ + matrixBuilder.newRowGroup(currentRow); + matrixBuilder.addNextValue(currentRow, sinkState, storm::utility::one()); + choiceLabeling.emplace_back(labelAll); + ++currentRow; + } + + return std::pair, std::vector>>(matrixBuilder.build(), std::move(choiceLabeling)); + } + + template + std::pair,std::vector>> SparseDtmcRegionModelChecker::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const{ + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Instantiation of flexible matrix is not supported for this type"); + } + +#ifdef STORM_HAVE_CARL + + template<> + void SparseDtmcRegionModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities){ + + if(true){ // eliminate all states with constant outgoing transitions + storm::storage::BitVector statesToEliminate = ~initialstates; + //todo: ordering of states important? + boost::optional> missingStateRewards; + for (auto const& state : statesToEliminate) { + bool onlyConstantOutgoingTransitions=true; + for(auto const& entry : flexibleMatrix.getRow(state)){ + if(!entry.getValue().isConstant()){ + onlyConstantOutgoingTransitions=false; + break; + } + } + if(onlyConstantOutgoingTransitions){ + eliminationModelChecker.eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); + subsystem.set(state,false); + } + } + + //Note: we could also "eliminate" the initial state to get rid of its selfloop + } + else if(false){ //eliminate all states with standard state elimination + boost::optional> missingStateRewards; + + storm::storage::BitVector statesToEliminate = ~initialstates; + std::vector states(statesToEliminate.begin(), statesToEliminate.end()); + + if (statePriorities) { + std::sort(states.begin(), states.end(), [&statePriorities] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statePriorities.get()[a] < statePriorities.get()[b]; }); + } + + STORM_LOG_DEBUG("Eliminating " << states.size() << " states using the state elimination technique." << std::endl); + for (auto const& state : states) { + eliminationModelChecker.eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); + } + subsystem=~statesToEliminate; + + } + else if(false){ //hybrid method + boost::optional> missingStateRewards; + storm::storage::BitVector statesToEliminate = ~initialstates; + 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(flexibleMatrix, oneStepProbabilities, initialstates, statesToEliminate, forwardTransitions, flexibleBackwardTransitions, 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(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions, missingStateRewards); + } + } + subsystem=~statesToEliminate; + } + std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << "states." << std::endl; + STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " states." << std::endl); + } + + template + void SparseDtmcRegionModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities){ + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "elimination of states not suported for this type"); + } + + + template<> + void SparseDtmcRegionModelChecker::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities){ + carl::VariablePool& varPool = carl::VariablePool::getInstance(); + + //first add a state variable for every state in the subsystem, providing that such a variable does not already exist. + for (storm::storage::sparse::state_type state : subsystem){ + if(stateProbVars[state].isZero()){ //variable does not exist yet + storm::Variable stateVar = varPool.getFreshVariable("p_" + std::to_string(state)); + std::shared_ptr>> cache(new carl::Cache>()); + storm::RationalFunction::PolyType stateVarAsPoly(storm::RationalFunction::PolyType::PolyType(stateVar), cache); + + //each variable is in the interval [0,1] + solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::GEQ, storm::RationalFunction(0)); + solver.add(storm::RationalFunction(stateVarAsPoly), storm::CompareRelation::LEQ, storm::RationalFunction(1)); + stateProbVars[state] = stateVarAsPoly; + } + } + + //now lets add the actual transitions + for (storm::storage::sparse::state_type state : subsystem){ + storm::RationalFunction reachProbability(oneStepProbabilities[state]); + for(auto const& transition : flexibleMatrix.getRow(state)){ + reachProbability += transition.getValue() * stateProbVars[transition.getColumn()]; + } + //Todo: depending on the objective (i.e. the formlua) it suffices to use LEQ or GEQ here... maybe this is faster? + solver.add(storm::RationalFunction(stateProbVars[state]), storm::CompareRelation::EQ, reachProbability); + } + } + + template + void SparseDtmcRegionModelChecker::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities){ + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "SMT formulation is not supported for this type"); + } + + template<> + void SparseDtmcRegionModelChecker::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compType){ + //We are going to build a new (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter + + //todo invent something better to obtain the substitutions. + //this only works as long as there is only one parameter per state, + // also: check whether the terms are linear/monotone(?) + + STORM_LOG_WARN("the probability restriction is experimental.. only works on linear terms with one parameter per state"); + storm::storage::sparse::state_type const numOfStates=subsystem.getNumberOfSetBits() + 2; //subsystem + target state + sink state + storm::models::sparse::StateLabeling stateLabeling(numOfStates); + stateLabeling.addLabel("init", storm::storage::BitVector(numOfStates, true)); + storm::storage::BitVector targetLabel(numOfStates, false); + targetLabel.set(numOfStates-2, true); + stateLabeling.addLabel("target", std::move(targetLabel)); + storm::storage::BitVector sinkLabel(numOfStates, false); + sinkLabel.set(numOfStates-1, true); + stateLabeling.addLabel("sink", std::move(sinkLabel)); + + std::map substitutionLB; + for(auto const& parRegion : regions){ + substitutionLB.insert(std::pair(parRegion.variable, parRegion.lowerBound)); + } + std::map substitutionUB; + for(auto const& parRegion : regions){ + substitutionUB.insert(std::pair(parRegion.variable, parRegion.upperBound)); + } + std::vector> substitutions; + substitutions.push_back(substitutionLB); + substitutions.push_back(substitutionUB); + + std::pair,std::vector>> instantiation = instantiateFlexibleMatrix(flexibleMatrix, substitutions, subsystem, true, oneStepProbabilities, true); + boost::optional> noStateRewards; + boost::optional> noTransitionRewards; + storm::models::sparse::Mdp mdp(instantiation.first, std::move(stateLabeling),noStateRewards,noTransitionRewards,instantiation.second); + + //we need the correct optimalityType for model checking as well as the correct relation for smt solving + storm::logic::OptimalityType opType; + storm::CompareRelation boundRelation; + switch (compType){ + case storm::logic::ComparisonType::Greater: + opType=storm::logic::OptimalityType::Minimize; + boundRelation=storm::CompareRelation::GEQ; + break; + case storm::logic::ComparisonType::GreaterEqual: + opType=storm::logic::OptimalityType::Minimize; + boundRelation=storm::CompareRelation::GEQ; + break; + case storm::logic::ComparisonType::Less: + opType=storm::logic::OptimalityType::Maximize; + boundRelation=storm::CompareRelation::LEQ; + break; + case storm::logic::ComparisonType::LessEqual: + opType=storm::logic::OptimalityType::Maximize; + boundRelation=storm::CompareRelation::LEQ; + break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); + } + + //perform model checking on the mdp + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); + storm::modelchecker::SparseMdpPrctlModelChecker modelChecker(mdp); + std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType); + std::vector resultVector = resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); + + //todo this is experimental.. + if (true){ + std::cout << "the matrix has " << mdp.getTransitionMatrix().getRowGroupCount() << "row groups, " << mdp.getTransitionMatrix().getRowCount() << " rows, and " << mdp.getTransitionMatrix().getColumnCount() << "columns" << std::endl; + boost::container::flat_set< uint_fast64_t> lbChoiceLabel; + lbChoiceLabel.insert(1); + lbChoiceLabel.insert(2); + storm::models::sparse::Dtmc dtmc(mdp.restrictChoiceLabels(lbChoiceLabel).getTransitionMatrix(), std::move(stateLabeling)); + //modelchecking on dtmc + storm::modelchecker::SparseDtmcPrctlModelChecker dtmcModelChecker(dtmc); + std::unique_ptr resultPtrDtmc = dtmcModelChecker.computeEventuallyProbabilities(eventuallyFormula); + std::vector resultVectorDtmc = resultPtrDtmc->asExplicitQuantitativeCheckResult().getValueVector(); + std::cout << "dtmc result with lower bounds:" << resultVectorDtmc[0]; + std::cout << "mdp result:" << resultVector[0]; + } + + + //formulate constraints for the solver + uint_fast64_t boundDenominator = 1.0/storm::settings::generalSettings().getPrecision(); //we need to approx. the obtained bounds as rational numbers + storm::storage::sparse::state_type subsystemState=0; //the subsystem uses other state indices + for(storm::storage::sparse::state_type state : subsystem){ + uint_fast64_t boundNumerator = resultVector[subsystemState]*boundDenominator; + storm::RationalFunction bound(boundNumerator); + bound = bound/boundDenominator; + //Todo: non-exact values might be problematic here... + solver.add(storm::RationalFunction(stateProbVars[state]), boundRelation, bound); + ++subsystemState; + } + } + + template + void SparseDtmcRegionModelChecker::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compType){ + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "restricting Probability Variables is not supported for this type"); + } + + template<> + bool SparseDtmcRegionModelChecker::checkRegion(storm::logic::Formula const& formula, std::vector::ParameterRegion> parameterRegions){ + //Note: this is an 'experimental' implementation + + std::chrono::high_resolution_clock::time_point timeStart = std::chrono::high_resolution_clock::now(); + + //Start with some preprocessing (inspired by computeUntilProbabilities...) + //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ]) + //get the (sub)formulae and the vector of target states + STORM_LOG_THROW(formula.isStateFormula(), storm::exceptions::IllegalArgumentException, "expected a stateFormula"); + STORM_LOG_THROW(formula.asStateFormula().isProbabilityOperatorFormula(), storm::exceptions::IllegalArgumentException, "expected a probabilityOperatorFormula"); + storm::logic::ProbabilityOperatorFormula const& probOpForm=formula.asStateFormula().asProbabilityOperatorFormula(); + STORM_LOG_THROW(probOpForm.hasBound(), storm::exceptions::IllegalArgumentException, "The formula has no bound"); + STORM_LOG_THROW(probOpForm.getSubformula().asPathFormula().isEventuallyFormula(), storm::exceptions::IllegalArgumentException, "expected an eventually subformula"); + storm::logic::EventuallyFormula const& eventuallyFormula = probOpForm.getSubformula().asPathFormula().asEventuallyFormula(); + std::unique_ptr targetStatesResultPtr = eliminationModelChecker.check(eventuallyFormula.getSubformula()); + storm::storage::BitVector const& targetStates = targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector(); + // Do some sanity checks to establish some required properties. + STORM_LOG_THROW(model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); + storm::storage::sparse::state_type initialState = *model.getInitialStates().begin(); + // 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 return the result. + if (model.getInitialStates().isDisjointFrom(maybeStates)) { + STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step."); + double res= statesWithProbability0.get(*model.getInitialStates().begin()) ? 0.0 : 1.0; + switch (probOpForm.getComparisonType()){ + case storm::logic::ComparisonType::Greater: + return (res > probOpForm.getBound()); + case storm::logic::ComparisonType::GreaterEqual: + return (res >= probOpForm.getBound()); + case storm::logic::ComparisonType::Less: + return (res < probOpForm.getBound()); + case storm::logic::ComparisonType::LessEqual: + return (res <= probOpForm.getBound()); + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula 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(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. + std::vector oneStepProbabilities = model.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); + // Determine the set of initial states of the sub-model. + storm::storage::BitVector newInitialStates = model.getInitialStates() % maybeStates; + // 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(); + + std::vector statePriorities = eliminationModelChecker.getStatePriorities(submatrix, submatrixTransposed, newInitialStates, oneStepProbabilities); + + // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. + FlexibleMatrix flexibleMatrix = eliminationModelChecker.getFlexibleSparseMatrix(submatrix); + FlexibleMatrix flexibleBackwardTransitions = eliminationModelChecker.getFlexibleSparseMatrix(submatrixTransposed, true); + + std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now(); + + // 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); + eliminateStates(subsystem, flexibleMatrix, oneStepProbabilities, flexibleBackwardTransitions, newInitialStates, submatrix, statePriorities); + + std::chrono::high_resolution_clock::time_point timeStateElemEnd = std::chrono::high_resolution_clock::now(); + + // SMT formulation of resulting pdtmc + storm::expressions::ExpressionManager manager; //this manager will do nothing as we will use carl expressions + storm::solver::Smt2SmtSolver solver(manager, true); + // we will introduce a variable for every state which encodes the probability to reach a target state from this state. + // we will store them as polynomials to easily use operations with rational functions + std::vector stateProbVars(subsystem.size(), storm::RationalFunction::PolyType(0)); + // todo maybe introduce the parameters already at this point? + formulateModelWithSMT(solver, stateProbVars, subsystem, flexibleMatrix, oneStepProbabilities); + + //the property should be satisfied in the initial state for all parameters. + //this is equivalent to: + //the negation of the property should not be satisfied for some parameter valuation. + //Hence, we flip the comparison relation and later check whether all the constraints are unsat. + storm::CompareRelation propertyCompRel; + switch (probOpForm.getComparisonType()){ + case storm::logic::ComparisonType::Greater: + propertyCompRel=storm::CompareRelation::LEQ; + break; + case storm::logic::ComparisonType::GreaterEqual: + propertyCompRel=storm::CompareRelation::LT; + break; + case storm::logic::ComparisonType::Less: + propertyCompRel=storm::CompareRelation::GEQ; + break; + case storm::logic::ComparisonType::LessEqual: + propertyCompRel=storm::CompareRelation::GT; + break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); + } + uint_fast64_t thresholdDenominator = 1.0/storm::settings::generalSettings().getPrecision(); + uint_fast64_t thresholdNumerator = probOpForm.getBound()*thresholdDenominator; + storm::RationalFunction threshold(thresholdNumerator); + threshold = threshold / thresholdDenominator; + solver.add(storm::RationalFunction(stateProbVars[*newInitialStates.begin()]), propertyCompRel, threshold); + + //the bounds for the parameters + solver.push(); + for(auto param : parameterRegions){ + storm::RawPolynomial lB(param.variable); + lB -= param.lowerBound; + solver.add(carl::Constraint(lB,storm::CompareRelation::GEQ)); + storm::RawPolynomial uB(param.variable); + uB -= param.upperBound; + solver.add(carl::Constraint(uB,storm::CompareRelation::LEQ)); + } + + std::chrono::high_resolution_clock::time_point timeSmtFormulationEnd = std::chrono::high_resolution_clock::now(); + + // find further restriction on probabilities + restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions, storm::logic::ComparisonType::Less); //probOpForm.getComparisonType()); + restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions, storm::logic::ComparisonType::Greater); + + std::chrono::high_resolution_clock::time_point timeRestrictingEnd = std::chrono::high_resolution_clock::now(); + + std::cout << "start solving ..." << std::endl; + bool result; + switch (solver.check()){ + case storm::solver::SmtSolver::CheckResult::Sat: + std::cout << "sat!" << std::endl; + result=false; + break; + case storm::solver::SmtSolver::CheckResult::Unsat: + std::cout << "unsat!" << std::endl; + result=true; + break; + case storm::solver::SmtSolver::CheckResult::Unknown: + std::cout << "unknown!" << std::endl; + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (Check-result: Unknown)") + result=false; + break; + default: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Could not solve the SMT-Problem (unexpected check-result)") + result=false; + } + + std::chrono::high_resolution_clock::time_point timeSolvingEnd = std::chrono::high_resolution_clock::now(); + + std::chrono::high_resolution_clock::duration timePreprocessing = timePreprocessingEnd - timeStart; + std::chrono::high_resolution_clock::duration timeStateElem = timeStateElemEnd - timePreprocessingEnd; + std::chrono::high_resolution_clock::duration timeSmtFormulation = timeSmtFormulationEnd - timeStateElemEnd; + std::chrono::high_resolution_clock::duration timeRestricting = timeRestrictingEnd - timeSmtFormulationEnd; + std::chrono::high_resolution_clock::duration timeSolving = timeSolvingEnd- timeRestrictingEnd; + std::chrono::high_resolution_clock::duration timeOverall = timeSolvingEnd - timeStart; + std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast(timePreprocessing); + std::chrono::milliseconds timeStateElemInMilliseconds = std::chrono::duration_cast(timeStateElem); + std::chrono::milliseconds timeSmtFormulationInMilliseconds = std::chrono::duration_cast(timeSmtFormulation); + std::chrono::milliseconds timeRestrictingInMilliseconds = std::chrono::duration_cast(timeRestricting); + std::chrono::milliseconds timeSolvingInMilliseconds = std::chrono::duration_cast(timeSolving); + std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast(timeOverall); + STORM_PRINT_AND_LOG(std::endl << "required time: " << timeOverallInMilliseconds.count() << "ms. Time Breakdown:" << std::endl); + STORM_PRINT_AND_LOG(" * " << timePreprocessingInMilliseconds.count() << "ms for Preprocessing" << std::endl); + STORM_PRINT_AND_LOG(" * " << timeStateElemInMilliseconds.count() << "ms for StateElemination" << std::endl); + STORM_PRINT_AND_LOG(" * " << timeSmtFormulationInMilliseconds.count() << "ms for SmtFormulation" << std::endl); + STORM_PRINT_AND_LOG(" * " << timeRestrictingInMilliseconds.count() << "ms for Restricting" << std::endl); + STORM_PRINT_AND_LOG(" * " << timeSolvingInMilliseconds.count() << "ms for Solving" << std::endl); + + return result; + } + + template + bool SparseDtmcRegionModelChecker::checkRegion(storm::logic::Formula const& formula, std::vector::ParameterRegion> parameterRegions){ + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Region check is not supported for this type"); + } + +#endif + +#ifdef STORM_HAVE_CARL + template class SparseDtmcRegionModelChecker; +#endif + } // namespace modelchecker +} // namespace storm diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h new file mode 100644 index 000000000..96f6385b2 --- /dev/null +++ b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h @@ -0,0 +1,89 @@ +#ifndef STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_ +#define STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_ + +//TODO remove useless includes + +#include "src/storage/sparse/StateType.h" +#include "src/models/sparse/Dtmc.h" +#include "src/models/sparse/Mdp.h" +//#include "src/modelchecker/AbstractModelChecker.h" +#include "src/utility/constants.h" +//#include "src/solver/SmtSolver.h" +#include "src/solver/Smt2SmtSolver.h" +#include "SparseDtmcEliminationModelChecker.h" + +namespace storm { + namespace modelchecker { + + template + class SparseDtmcRegionModelChecker { + public: + explicit SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc const& model); + + virtual bool canHandle(storm::logic::Formula const& formula) const; + +#ifdef STORM_HAVE_CARL + struct ParameterRegion{ + storm::Variable variable; + storm::RationalFunction::CoeffType lowerBound; + storm::RationalFunction::CoeffType upperBound; + }; + + + /*! + * Checks whether the given formula holds for all possible parameters that satisfy the given parameter regions + * ParameterRegions should contain all parameters (not mentioned parameters are assumed to be arbitrary reals) + */ + bool checkRegion(storm::logic::Formula const& formula, std::vector parameterRegions); +#endif + private: + typedef typename storm::modelchecker::SparseDtmcEliminationModelChecker::FlexibleSparseMatrix FlexibleMatrix; + + #ifdef STORM_HAVE_CARL + /*! + * Instantiates the matrix, i.e., evaluate the occurring functions according to the given substitutions of the variables. + * One row of the given flexible matrix will be one rowgroup in the returned matrix, consisting of one row for every substitution. + * The returned matrix can be seen as the transition matrix of an MDP with the action labeling given by the returned vector of sets. + * Only the rows selected by the given filter are considered. (filter should have size==matrix.getNumberOfRows()) + * An exception is thrown if there is a transition from a selected state to an unselected state + * If one step probabilities are given, a new state is added which can be considered as target state. + * The "missing" probability can be redirected to a sink state + * By convention, the target state will have index filter.getNumberOfSetBits() and the sink state will be the state with the highest index (so right after the target state) + * + * @param matrix the considered flexible matrix + * @param substitutions A list of mappings, each assigning a constant value to every variable + * @param filter selects the rows of this flexibleMatrix, that will be considered + * @param addSinkState adds a state with a self loop to which the "missing" probability will lead + * @param oneStepProbabilities if given, a new state is added to which there are transitions for all non-zero entries in this vector + * @param addSelfLoops if set, zero valued selfloops will be added in every row + * + * @return A matrix with constant (double) entries and a choice labeling + */ + std::pair,std::vector>> instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState=true, std::vector const& oneStepProbabilities=std::vector(), bool addSelfLoops=true) const; + //todo add const keyword + + //eliminates some of the states according to different strategies. + void eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities = {}); + + void formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities); + + void restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, std::vector const& regions, storm::logic::ComparisonType const& compTypeOfProperty); +#endif + + + + // The model this model checker is supposed to analyze. + storm::models::sparse::Dtmc const& model; + + // Instance of an elimination model checker to access its functions + storm::modelchecker::SparseDtmcEliminationModelChecker eliminationModelChecker; + + // A comparator that can be used to compare constants. + storm::utility::ConstantsComparator comparator; + + }; + + } // namespace modelchecker +} // namespace storm + +#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEDTMCREGIONMODELCHECKER_H_ */ diff --git a/src/utility/cli.h b/src/utility/cli.h index 5a71a5aa1..59d20e7eb 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -65,6 +65,7 @@ log4cplus::Logger printer; // Headers for model checking. #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h" +#include "src/modelchecker/reachability/SparseDtmcRegionModelChecker.h" #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h" @@ -475,21 +476,21 @@ namespace storm { if(settings.isParametricRegionSet()){ std::cout << std::endl; //experimental implementation! check some hardcoded region - std::vector::ParameterRegion> regions; + std::vector::ParameterRegion> regions; storm::RationalFunction::CoeffType zeroPointOne(1); zeroPointOne = zeroPointOne/10; - storm::modelchecker::SparseDtmcEliminationModelChecker::ParameterRegion param1; + storm::modelchecker::SparseDtmcRegionModelChecker::ParameterRegion param1; param1.lowerBound= zeroPointOne*zeroPointOne*78; param1.upperBound= zeroPointOne*zeroPointOne*82; param1.variable=carl::VariablePool::getInstance().findVariableWithName("pL"); regions.push_back(param1); - storm::modelchecker::SparseDtmcEliminationModelChecker::ParameterRegion param2; + storm::modelchecker::SparseDtmcRegionModelChecker::ParameterRegion param2; param2.lowerBound= zeroPointOne*zeroPointOne*78; param2.upperBound= zeroPointOne*zeroPointOne*82;; param2.variable=carl::VariablePool::getInstance().findVariableWithName("pK"); regions.push_back(param2); - storm::modelchecker::SparseDtmcEliminationModelChecker modelchecker(*dtmc); + storm::modelchecker::SparseDtmcRegionModelChecker modelchecker(*dtmc); bool result = modelchecker.checkRegion(*formula.get(), regions); std::cout << "... done." << std::endl; if (result){