diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index 0025c82a5..8d8a91302 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -17,6 +17,7 @@ #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/InvalidStateException.h" #include "exceptions/UnexpectedException.h" +#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h" namespace storm { namespace modelchecker { @@ -963,7 +964,7 @@ namespace storm { } template - bool SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::hasSelfLoop(storm::storage::sparse::state_type state) { + bool SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::hasSelfLoop(storm::storage::sparse::state_type state) const { for (auto const& entry : this->getRow(state)) { if (entry.getColumn() < state) { continue; @@ -989,32 +990,96 @@ namespace storm { #ifdef STORM_HAVE_CARL template<> - storm::storage::SparseMatrix SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::map substitutions){ + storm::storage::SparseMatrix SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::map const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const { - //Check if the CoeffType is as expected + //Check if the arguments are as expected STORM_LOG_THROW((std::is_same::value), storm::exceptions::IllegalArgumentException, "Unexpected Type of Coefficients"); - - //get a Matrix builder - index_type numElements=0; - for(row_type const& row : this->data){ - numElements += row.size(); + 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 data for a Matrix builder as well as a mapping from old state indices to the new ones + index_type numTransitions=0; + 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){ + numTransitions += this->getRow(state).size(); + if(addSelfLoops && !hasSelfLoop(state)){ + ++numTransitions; + } + if(!oneStepProbabilities.empty() && !oneStepProbabilities[state].isZero()){ + ++numTransitions; + } + if(addSinkState){ + ++numTransitions; //we always add a transition here.. Todo: consider other ways with less memory consumption to handle this + } + newStateIndexMap[state]=newStateIndex; + ++newStateIndex; } - storm::storage::SparseMatrixBuilder matrixBuilder(this->getNumberOfRows(), this->getNumberOfRows(), numElements); - - //fill in the data... - for(index_type rowIndex=0; rowIndex < this->getNumberOfRows(); ++rowIndex){ - for(auto const& entry : this->getRow(rowIndex)){ - double value = cln::double_approx(entry.getValue().evaluate(substitutions)); - matrixBuilder.addNextValue(rowIndex, entry.getColumn(), value); + 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; + ++numTransitions; + } + if(addSinkState){ + sinkState=numStates; + ++numStates; + ++numTransitions; + } + storm::storage::SparseMatrixBuilder matrixBuilder(numStates, numStates, numTransitions); + //fill in the data row by row + for(auto const& oldStateIndex : filter){ + double missingProbability = 1.0; + if(this->getRow(oldStateIndex).empty()){ + if(addSelfLoops){ + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], newStateIndexMap[oldStateIndex], 0.0); + } + } + else{ + const_iterator entry = this->getRow(oldStateIndex).begin(); + for(; entrygetRow(oldStateIndex).end() && entry->getColumn()getValue().evaluate(substitutions)); + missingProbability-=value; + storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; + STORM_LOG_THROW(columngetColumn()!=oldStateIndex){ + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], newStateIndexMap[oldStateIndex], 0.0); + } + for(; entry < this->getRow(oldStateIndex).end(); ++entry){ + double value = cln::double_approx(entry->getValue().evaluate(substitutions)); + missingProbability-=value; + storm::storage::sparse::state_type column = newStateIndexMap[entry->getColumn()]; + STORM_LOG_THROW(column.isZero(missingProbability)){ + STORM_LOG_ASSERT(missingProbability> -storm::settings::generalSettings().getPrecision(), "The missing probability is negative."); + matrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, missingProbability); } } + if (!oneStepProbabilities.empty()){ + matrixBuilder.addNextValue(targetState, targetState, 1.0); + } + if (addSinkState){ + matrixBuilder.addNextValue(sinkState, sinkState, 1.0); + } return matrixBuilder.build(); } template - storm::storage::SparseMatrix SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::map substitutions){ - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Instantiation of flexible matrix is not supported for this type"); + storm::storage::SparseMatrix SparseDtmcEliminationModelChecker::FlexibleSparseMatrix::instantiateAsDouble(std::map 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 @@ -1051,7 +1116,7 @@ namespace storm { template void SparseDtmcEliminationModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleSparseMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleSparseMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates){ //.. simple strategy for now... do not eliminate anything - /* + /* ... or all? boost::optional> missingStateRewards; storm::storage::BitVector statesToEliminate = ~initialstates; @@ -1101,6 +1166,67 @@ namespace storm { 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) DTMC + STORM_LOG_WARN("the probability restriction is not really correct, it only helps if there is a 'sat' answer"); + 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 substitutions; + for(auto const& parRegion : regions){ + substitutions.insert(std::pair(parRegion.variable, parRegion.upperBound)); //todo: (upper+lower)/2 ? + } + storm::models::sparse::Dtmc dtmc(flexibleMatrix.instantiateAsDouble(substitutions, subsystem, true, oneStepProbabilities, true), std::move(stateLabeling)); + + //perform model checking on this dtmc + storm::modelchecker::SparseDtmcPrctlModelChecker modelChecker(dtmc); + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); + std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula); + std::vector resultVector = resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); + + //formulate constraints for the solver + storm::CompareRelation boundRelation; + switch (compType){ + case storm::logic::ComparisonType::Greater: + boundRelation=storm::CompareRelation::LEQ; + break; + case storm::logic::ComparisonType::GreaterEqual: + boundRelation=storm::CompareRelation::LEQ; + break; + case storm::logic::ComparisonType::Less: + boundRelation=storm::CompareRelation::GEQ; + break; + case storm::logic::ComparisonType::LessEqual: + boundRelation=storm::CompareRelation::GEQ; + break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); + } + 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 @@ -1115,8 +1241,8 @@ namespace storm { 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& eventFormula = probOpForm.getSubformula().asPathFormula().asEventuallyFormula(); - std::unique_ptr targetStatesResultPtr = this->check(eventFormula.getSubformula()); + 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."); @@ -1215,7 +1341,8 @@ namespace storm { std::chrono::high_resolution_clock::time_point timeSmtFormulationEnd = std::chrono::high_resolution_clock::now(); - // TODO find further restriction on probabilities + // find further restriction on probabilities + restrictProbabilityVariables(solver,stateProbVars,subsystem,flexibleMatrix,oneStepProbabilities, parameterRegions, probOpForm.getComparisonType()); std::chrono::high_resolution_clock::time_point timeRestrictingEnd = std::chrono::high_resolution_clock::now(); diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h index ab595dcea..064a70499 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h @@ -68,17 +68,28 @@ namespace storm { * @param matrix The matrix in which to look for the loop. * @return True iff the given state has a self-loop with an arbitrary probability in the given probability matrix. */ - bool hasSelfLoop(storm::storage::sparse::state_type state); + 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 substitution of the variables + * Instantiates the matrix, i.e., evaluate the occurring functions according to the given substitution of the variables. + * 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 mapping that assigns 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 */ - storm::storage::SparseMatrix instantiateAsDouble(std::map substitutions); + storm::storage::SparseMatrix instantiateAsDouble(std::map 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: @@ -100,6 +111,8 @@ namespace storm { 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;