diff --git a/examples/pdtmc/crowds/crowds_10-5.pm b/examples/pdtmc/crowds/crowds_10-5.pm index dcf6e7930..ab5dbd438 100644 --- a/examples/pdtmc/crowds/crowds_10-5.pm +++ b/examples/pdtmc/crowds/crowds_10-5.pm @@ -10,7 +10,7 @@ dtmc // Model parameters const double PF; // forwarding probability -const double badC = 0.167; // probability that member is untrustworthy +const double badC; // probability that member is untrustworthy // Probability of forwarding // const double PF = 0.8; diff --git a/examples/pdtmc/die/die.pm b/examples/pdtmc/die/die.pm new file mode 100644 index 000000000..84e466a18 --- /dev/null +++ b/examples/pdtmc/die/die.pm @@ -0,0 +1,33 @@ +// Knuth's model of a fair die using only fair coins +dtmc + +const double coinProb; + +module die + + // local state + s : [0..7] init 0; + // value of the dice + d : [0..6] init 0; + + [] s=0 -> coinProb : (s'=1) + (1-coinProb) : (s'=2); + [] s=1 -> coinProb : (s'=3) + (1-coinProb) : (s'=4); + [] s=2 -> coinProb : (s'=5) + (1-coinProb) : (s'=6); + [] s=3 -> coinProb : (s'=1) + (1-coinProb) : (s'=7) & (d'=1); + [] s=4 -> coinProb : (s'=7) & (d'=2) + (1-coinProb) : (s'=7) & (d'=3); + [] s=5 -> coinProb : (s'=7) & (d'=4) + (1-coinProb) : (s'=7) & (d'=5); + [] s=6 -> coinProb : (s'=2) + (1-coinProb) : (s'=7) & (d'=6); + [] s=7 -> 1: (s'=7); + +endmodule + +rewards "coin_flips" + [] s<7 : 1; +endrewards + +label "one" = s=7&d=1; +//label "two" = s=7&d=2; +//label "three" = s=7&d=3; +//label "four" = s=7&d=4; +//label "five" = s=7&d=5; +//label "six" = s=7&d=6; diff --git a/src/adapters/ExplicitModelAdapter.h b/src/adapters/ExplicitModelAdapter.h index 35a7d1249..89dce2caf 100644 --- a/src/adapters/ExplicitModelAdapter.h +++ b/src/adapters/ExplicitModelAdapter.h @@ -29,6 +29,7 @@ #include "src/storage/SparseMatrix.h" #include "src/settings/SettingsManager.h" #include "src/utility/macros.h" +#include "src/utility/ConstantsComparator.h" #include "src/exceptions/WrongFormatException.h" #include "src/storage/expressions/ExpressionEvaluation.h" @@ -269,7 +270,7 @@ namespace storm { return result; } - static std::list> getUnlabeledTransitions(storm::prism::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue& stateQueue, expressions::ExpressionEvaluation& eval) { + static std::list> getUnlabeledTransitions(storm::prism::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue& stateQueue, storm::utility::ConstantsComparator const& comparator, expressions::ExpressionEvaluation& eval) { std::list> result; StateType const* currentState = stateInformation.reachableStates[stateIndex]; @@ -316,14 +317,14 @@ namespace storm { } // Check that the resulting distribution is in fact a distribution. - STORM_LOG_THROW(storm::utility::isOne(probabilitySum), storm::exceptions::WrongFormatException, "Probabilities do not sum to one for command '" << command << "' (sum = " << probabilitySum << ")."); + STORM_LOG_THROW(!comparator.isConstant(probabilitySum) || comparator.isOne(probabilitySum), storm::exceptions::WrongFormatException, "Probabilities do not sum to one for command '" << command << "' (sum = " << probabilitySum << ")."); } } return result; } - static std::list> getLabeledTransitions(storm::prism::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue& stateQueue, expressions::ExpressionEvaluation& eval) { + static std::list> getLabeledTransitions(storm::prism::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue& stateQueue, storm::utility::ConstantsComparator const& comparator, expressions::ExpressionEvaluation& eval) { std::list> result; for (std::string const& action : program.getActions()) { @@ -422,10 +423,7 @@ namespace storm { } // Check that the resulting distribution is in fact a distribution. - if (!storm::utility::isOne(probabilitySum)) { - LOG4CPLUS_ERROR(logger, "Sum of update probabilities do not some to one for some command."); - throw storm::exceptions::WrongFormatException() << "Sum of update probabilities do sum to " << probabilitySum << " to one for some command."; - } + STORM_LOG_THROW(!comparator.isConstant(probabilitySum) || comparator.isOne(probabilitySum), storm::exceptions::WrongFormatException, "Sum of update probabilities do sum to " << probabilitySum << " to one for some command."); // Dispose of the temporary maps. delete currentTargetStates; @@ -467,7 +465,7 @@ namespace storm { * @return A tuple containing a vector with all rows at which the nondeterministic choices of each state begin * and a vector containing the labels associated with each choice. */ - static std::vector> buildMatrices(storm::prism::Program const& program, VariableInformation const& variableInformation, std::vector const& transitionRewards, StateInformation& stateInformation, bool deterministicModel, storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, storm::storage::SparseMatrixBuilder& transitionRewardMatrixBuilder, expressions::ExpressionEvaluation& eval) { + static std::vector> buildMatrices(storm::prism::Program const& program, VariableInformation const& variableInformation, std::vector const& transitionRewards, StateInformation& stateInformation, bool deterministicModel, storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, storm::storage::SparseMatrixBuilder& transitionRewardMatrixBuilder, storm::utility::ConstantsComparator const& comparator, expressions::ExpressionEvaluation& eval) { std::vector> choiceLabels; // Initialize a queue and insert the initial state. @@ -498,8 +496,8 @@ namespace storm { uint_fast64_t currentState = stateQueue.front(); // Retrieve all choices for the current state. - std::list> allUnlabeledChoices = getUnlabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue, eval); - std::list> allLabeledChoices = getLabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue, eval); + std::list> allUnlabeledChoices = getUnlabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue, comparator, eval); + std::list> allLabeledChoices = getLabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue, comparator, eval); uint_fast64_t totalNumberOfChoices = allUnlabeledChoices.size() + allLabeledChoices.size(); @@ -646,7 +644,6 @@ namespace storm { ModelComponents modelComponents; expressions::ExpressionEvaluation eval; - VariableInformation variableInformation; for (auto const& integerVariable : program.getGlobalIntegerVariables()) { variableInformation.variableToBoundsMap[integerVariable.getName()] = std::make_pair(integerVariable.getLowerBoundExpression().evaluateAsInt(), integerVariable.getUpperBoundExpression().evaluateAsInt()); @@ -670,7 +667,8 @@ namespace storm { // Build the transition and reward matrices. storm::storage::SparseMatrixBuilder transitionMatrixBuilder(0, 0, 0, false, !deterministicModel, 0); storm::storage::SparseMatrixBuilder transitionRewardMatrixBuilder(0, 0, 0, false, !deterministicModel, 0); - modelComponents.choiceLabeling = buildMatrices(program, variableInformation, rewardModel.getTransitionRewards(), stateInformation, deterministicModel, transitionMatrixBuilder, transitionRewardMatrixBuilder, eval); + storm::utility::ConstantsComparator comparator; + modelComponents.choiceLabeling = buildMatrices(program, variableInformation, rewardModel.getTransitionRewards(), stateInformation, deterministicModel, transitionMatrixBuilder, transitionRewardMatrixBuilder, comparator, eval); // Finalize the resulting matrices. modelComponents.transitionMatrix = transitionMatrixBuilder.build(); diff --git a/src/adapters/extendedCarl.h b/src/adapters/extendedCarl.h index 90c1cef4d..fddb4950f 100644 --- a/src/adapters/extendedCarl.h +++ b/src/adapters/extendedCarl.h @@ -1,4 +1,4 @@ -/** +/** * @file: extendedCarl.h * @author: Sebastian Junges * @@ -13,22 +13,29 @@ #include #include #include +#include namespace carl { -template -inline size_t hash_value(carl::MultivariatePolynomial const& p) -{ - std::hash> h; - return h(p); -} -template -inline size_t hash_value(carl::RationalFunction const& f) -{ - std::hash h; - return h(f.nominator()) ^ h(f.denominator()); -} - + template + inline size_t hash_value(carl::MultivariatePolynomial const& p) + { + std::hash> h; + return h(p); + } + template + inline size_t hash_value(carl::RationalFunction const& f) + { + std::hash h; + return h(f.nominator()) ^ h(f.denominator()); + } + + template + inline size_t hash_value(carl::FactorizedPolynomial const& p) + { + std::hash> h; + return h(p); + } } diff --git a/src/modelchecker/reachability/DirectEncoding.h b/src/modelchecker/reachability/DirectEncoding.h index 2c9f43e7c..69ddeeb28 100644 --- a/src/modelchecker/reachability/DirectEncoding.h +++ b/src/modelchecker/reachability/DirectEncoding.h @@ -1,4 +1,4 @@ -/** +/** * @file: DirectEncoding.h * @author: Sebastian Junges * @@ -12,96 +12,69 @@ namespace storm { - namespace modelchecker - { - namespace reachability - { - class DirectEncoding - { - public: - template - std::string encodeAsSmt2(const storm::models::Dtmc& model, std::vector parameters, storm::storage::BitVector initialStates, storm::storage::BitVector finalStates, const typename T::CoeffType& threshold, bool lessequal = false) - { - carl::io::WriteTosmt2Stream smt2; - uint_fast64_t nrStates = model.getNumberOfStates(); - carl::VariablePool& vpool = carl::VariablePool::getInstance(); - std::vector stateVars; - for(carl::Variable p : parameters) - { - smt2 << ("parameter_bound_" + vpool.getName(p)); - smt2 << carl::io::smt2node::AND; - smt2 << carl::Constraint(Polynomial(p), carl::CompareRelation::GT); - smt2 << carl::Constraint(Polynomial(p) - Polynomial(1), carl::CompareRelation::LT); - smt2 << carl::io::smt2node::CLOSENODE; - } - for(uint_fast64_t state = 0; state < nrStates-1; ++state) - { - - carl::Variable stateVar = vpool.getFreshVariable("s_" + std::to_string(state)); - stateVars.push_back(stateVar); - if(!finalStates[state]) - { - smt2 << ("state_bound_" + std::to_string(state)); - smt2 << carl::io::smt2node::AND; - smt2 << carl::Constraint(Polynomial(stateVar), carl::CompareRelation::GT); - smt2 << carl::Constraint(Polynomial(stateVar) - Polynomial(1), carl::CompareRelation::LT); - smt2 << carl::io::smt2node::CLOSENODE; - } - - } - - smt2.setAutomaticLineBreaks(true); - Polynomial initStateReachSum; - for(uint_fast64_t state = 0; state < nrStates-1; ++state) - { - if(initialStates[state]) - { - initStateReachSum += stateVars[state]; - } - if(finalStates[state]) - { - //smt2 << carl::Constraint(Polynomial(stateVars[state]) - Polynomial(1), carl::CompareRelation::EQ); - } - else - { - T reachpropPol(0); - for(auto const& transition : model.getRows(state)) - { - if(finalStates[transition.first]) - { - reachpropPol += transition.second; - } - else if(transition.first == nrStates - 1) - { - // intentionally empty. - } - else - { - reachpropPol += transition.second * stateVars[transition.first]; - } - - } - smt2 << ("transition_" + std::to_string(state)); - smt2 << carl::Constraint(reachpropPol - stateVars[state], carl::CompareRelation::EQ); - } - } - //smt2 << carl::Constraint(Polynomial(stateVars[nrStates-1]), carl::CompareRelation::EQ); - - - smt2 << ("reachability"); - - carl::CompareRelation thresholdRelation = lessequal ? carl::CompareRelation::LEQ : carl::CompareRelation::GEQ; - smt2 << carl::Constraint(initStateReachSum - threshold, thresholdRelation); - smt2 << carl::io::smt2flag::CHECKSAT; - smt2 << carl::io::smt2flag::MODEL; - smt2 << carl::io::smt2flag::UNSAT_CORE; - std::stringstream strm; - strm << smt2; - return strm.str(); - } - }; - } - } + namespace modelchecker + { + namespace reachability + { + class DirectEncoding + { + public: + template + std::string encodeAsSmt2(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& oneStepProbabilities, std::set const& parameters, storm::storage::BitVector const& initialStates, typename T::CoeffType const& threshold, bool lessequal = false) { + + carl::io::WriteTosmt2Stream smt2; + uint_fast64_t numberOfStates = transitionMatrix.getRowCount(); + carl::VariablePool& vpool = carl::VariablePool::getInstance(); + std::vector stateVars; + for (carl::Variable const& p : parameters) { + smt2 << ("parameter_bound_" + vpool.getName(p)); + smt2 << carl::io::smt2node::AND; + smt2 << carl::Constraint(Polynomial::PolyType(p), carl::CompareRelation::GT); + smt2 << carl::Constraint(Polynomial::PolyType(p) - Polynomial::PolyType(1), carl::CompareRelation::LT); + smt2 << carl::io::smt2node::CLOSENODE; + } + + for (uint_fast64_t state = 0; state < numberOfStates; ++state) { + carl::Variable stateVar = vpool.getFreshVariable("s_" + std::to_string(state)); + stateVars.push_back(stateVar); + smt2 << ("state_bound_" + std::to_string(state)); + smt2 << carl::io::smt2node::AND; + smt2 << carl::Constraint(Polynomial::PolyType(stateVar), carl::CompareRelation::GT); + smt2 << carl::Constraint(Polynomial::PolyType(stateVar) - Polynomial::PolyType(1), carl::CompareRelation::LT); + smt2 << carl::io::smt2node::CLOSENODE; + } + + smt2.setAutomaticLineBreaks(true); + Polynomial::PolyType initStateReachSum; + for (uint_fast64_t state = 0; state < numberOfStates; ++state) { + T reachpropPol; + for (auto const& transition : transitionMatrix.getRow(state)) { +// reachpropPol += transition.getValue() * stateVars[transition.getColumn()]; + } + reachpropPol += oneStepProbabilities[state]; + smt2 << ("transition_" + std::to_string(state)); +// smt2 << carl::Constraint(reachpropPol - stateVars[state], carl::CompareRelation::EQ); + } + + smt2 << ("reachability"); + + carl::CompareRelation thresholdRelation = lessequal ? carl::CompareRelation::LEQ : carl::CompareRelation::GEQ; + smt2 << carl::io::smt2node::OR; + for (uint_fast64_t state : initialStates) { + smt2 << carl::Constraint(Polynomial::PolyType(stateVars[state]) - threshold, thresholdRelation); + } + smt2 << carl::io::smt2node::CLOSENODE; + + smt2 << carl::io::smt2flag::CHECKSAT; + smt2 << carl::io::smt2flag::MODEL; + smt2 << carl::io::smt2flag::UNSAT_CORE; + std::stringstream strm; + strm << smt2; + return strm.str(); + } + }; + } + } } #endif \ No newline at end of file diff --git a/src/modelchecker/reachability/SparseSccModelChecker.cpp b/src/modelchecker/reachability/SparseSccModelChecker.cpp index 51ba38c64..7aae3f0b4 100644 --- a/src/modelchecker/reachability/SparseSccModelChecker.cpp +++ b/src/modelchecker/reachability/SparseSccModelChecker.cpp @@ -10,37 +10,64 @@ #include "src/utility/graph.h" #include "src/utility/vector.h" #include "src/utility/macros.h" +#include "src/utility/ConstantsComparator.h" namespace storm { namespace modelchecker { namespace reachability { template - static ValueType&& simplify(ValueType&& value) { - // In the general case, we don't to anything here, but merely return the value. If something else is - // supposed to happen here, the templated function can be specialized for this particular type. - return std::forward(value); - } - - static RationalFunction&& simplify(RationalFunction&& value) { - // In the general case, we don't to anything here, but merely return the value. If something else is - // supposed to happen here, the templated function can be specialized for this particular type. - STORM_LOG_TRACE("Simplifying " << value << " ... "); - value.simplify(); - STORM_LOG_TRACE("done."); - return std::forward(value); - } - - template - static storm::storage::MatrixEntry&& simplify(storm::storage::MatrixEntry&& matrixEntry) { - simplify(matrixEntry.getValue()); - return std::move(matrixEntry); - } - - template - static storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry) { - matrixEntry.setValue(simplify(matrixEntry.getValue())); - return matrixEntry; + ValueType SparseSccModelChecker::computeReachabilityProbability(storm::storage::SparseMatrix const& transitionMatrix, std::vector& oneStepProbabilities, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional> const& distances) { + auto totalTimeStart = std::chrono::high_resolution_clock::now(); + auto conversionStart = std::chrono::high_resolution_clock::now(); + + // Create a bit vector that represents the subsystem of states we still have to eliminate. + storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount(), true); + + // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily. + FlexibleSparseMatrix flexibleMatrix = getFlexibleSparseMatrix(transitionMatrix); + FlexibleSparseMatrix flexibleBackwardTransitions = getFlexibleSparseMatrix(backwardTransitions, true); + auto conversionEnd = std::chrono::high_resolution_clock::now(); + + // Then, we recursively treat all SCCs. + storm::utility::ConstantsComparator comparator; + auto modelCheckingStart = std::chrono::high_resolution_clock::now(); + std::vector entryStateQueue; + uint_fast64_t maximalDepth = treatScc(flexibleMatrix, oneStepProbabilities, initialStates, subsystem, transitionMatrix, flexibleBackwardTransitions, false, 0, storm::settings::parametricSettings().getMaximalSccSize(), entryStateQueue, comparator, distances); + + // 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::parametricSettings().isEliminateEntryStatesLastSet()) { + for (auto const& state : entryStateQueue) { + eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions); + } + } + + // Make sure that we have eliminated all transitions from the initial state. + STORM_LOG_ASSERT(flexibleMatrix.getRow(*initialStates.begin()).empty(), "The transitions of the initial states are non-empty."); + + auto modelCheckingEnd = std::chrono::high_resolution_clock::now(); + auto totalTimeEnd = std::chrono::high_resolution_clock::now(); + + if (storm::settings::generalSettings().isShowStatisticsSet()) { + auto conversionTime = conversionEnd - conversionStart; + auto modelCheckingTime = modelCheckingEnd - modelCheckingStart; + auto totalTime = totalTimeEnd - totalTimeStart; + + STORM_PRINT_AND_LOG(std::endl); + STORM_PRINT_AND_LOG("Time breakdown:" << std::endl); + STORM_PRINT_AND_LOG(" * time for conversion: " << std::chrono::duration_cast(conversionTime).count() << "ms" << std::endl); + STORM_PRINT_AND_LOG(" * time for checking: " << std::chrono::duration_cast(modelCheckingTime).count() << "ms" << std::endl); + STORM_PRINT_AND_LOG("------------------------------------------" << std::endl); + STORM_PRINT_AND_LOG(" * total time: " << std::chrono::duration_cast(totalTime).count() << "ms" << std::endl); + STORM_PRINT_AND_LOG(std::endl); + STORM_PRINT_AND_LOG("Other:" << std::endl); + STORM_PRINT_AND_LOG(" * number of states eliminated: " << transitionMatrix.getRowCount() << std::endl); + STORM_PRINT_AND_LOG(" * maximal depth of SCC decomposition: " << maximalDepth << std::endl); + } + + // Now, we return the value for the only initial state. + return storm::utility::simplify(oneStepProbabilities[*initialStates.begin()]); } template @@ -74,30 +101,26 @@ namespace storm { // Do some sanity checks to establish some required properties. STORM_LOG_THROW(dtmc.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); - storm::storage::sparse::state_type initialStateIndex = *dtmc.getInitialStates().begin(); // Then, compute the subset of states that has a probability of 0 or 1, respectively. - auto totalTimeStart = std::chrono::high_resolution_clock::now(); - auto preprocessingStart = std::chrono::high_resolution_clock::now(); std::pair statesWithProbability01 = storm::utility::graph::performProb01(dtmc, phiStates, psiStates); 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 (!maybeStates.get(initialStateIndex)) { - return statesWithProbability0.get(initialStateIndex) ? storm::utility::constantZero() : storm::utility::constantOne(); + if (dtmc.getInitialStates().isDisjointFrom(maybeStates)) { + STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step."); + return statesWithProbability0.get(*dtmc.getInitialStates().begin()) ? storm::utility::constantZero() : storm::utility::constantOne(); } - + // 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(dtmc.getTransitionMatrix(), dtmc.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; - auto preprocessingEnd = std::chrono::high_resolution_clock::now(); - + // Create a vector for the probabilities to go to a state with probability 1 in one step. - auto conversionStart = std::chrono::high_resolution_clock::now(); std::vector oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); // Determine the set of initial states of the sub-DTMC. @@ -107,17 +130,19 @@ namespace storm { storm::storage::SparseMatrix submatrix = dtmc.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates); // To be able to apply heuristics later, we now determine the distance of each state to the initial state. + // We start by setting up the data structures. std::vector> stateQueue; - stateQueue.reserve(maybeStates.getNumberOfSetBits()); - std::vector distances(maybeStates.getNumberOfSetBits()); - storm::storage::BitVector statesInQueue(maybeStates.getNumberOfSetBits()); + stateQueue.reserve(submatrix.getRowCount()); + storm::storage::BitVector statesInQueue(submatrix.getRowCount()); + std::vector distances(submatrix.getRowCount()); + storm::storage::sparse::state_type currentPosition = 0; for (auto const& initialState : newInitialStates) { stateQueue.emplace_back(initialState, 0); statesInQueue.set(initialState); } - // Perform a BFS. + // And then perform the BFS. while (currentPosition < stateQueue.size()) { std::pair const& stateDistancePair = stateQueue[currentPosition]; distances[stateDistancePair.first] = stateDistancePair.second; @@ -131,51 +156,7 @@ namespace storm { ++currentPosition; } - // Create a bit vector that represents the subsystem of states we still have to eliminate. - storm::storage::BitVector subsystem = storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true); - - // 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(submatrix.transpose(), true); - auto conversionEnd = std::chrono::high_resolution_clock::now(); - - // Then, we recursively treat all SCCs. - auto modelCheckingStart = std::chrono::high_resolution_clock::now(); - std::vector entryStateQueue; - uint_fast64_t maximalDepth = treatScc(dtmc, flexibleMatrix, oneStepProbabilities, newInitialStates, subsystem, submatrix, flexibleBackwardTransitions, false, 0, storm::settings::parametricSettings().getMaximalSccSize(), entryStateQueue, distances); - - // 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::parametricSettings().isEliminateEntryStatesLastSet()) { - for (auto const& state : entryStateQueue) { - eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions); - } - } - - auto modelCheckingEnd = std::chrono::high_resolution_clock::now(); - auto totalTimeEnd = std::chrono::high_resolution_clock::now(); - - if (storm::settings::generalSettings().isShowStatisticsSet()) { - auto preprocessingTime = preprocessingEnd - preprocessingStart; - auto conversionTime = conversionEnd - conversionStart; - auto modelCheckingTime = modelCheckingEnd - modelCheckingStart; - auto totalTime = totalTimeEnd - totalTimeStart; - - STORM_PRINT_AND_LOG(std::endl); - STORM_PRINT_AND_LOG("Time breakdown:" << std::endl); - STORM_PRINT_AND_LOG(" * time for preprocessing: " << std::chrono::duration_cast(preprocessingTime).count() << "ms" << std::endl); - STORM_PRINT_AND_LOG(" * time for conversion: " << std::chrono::duration_cast(conversionTime).count() << "ms" << std::endl); - STORM_PRINT_AND_LOG(" * time for checking: " << std::chrono::duration_cast(modelCheckingTime).count() << "ms" << std::endl); - STORM_PRINT_AND_LOG("------------------------------------------" << std::endl); - STORM_PRINT_AND_LOG(" * total time: " << std::chrono::duration_cast(totalTime).count() << "ms" << std::endl); - STORM_PRINT_AND_LOG(std::endl); - STORM_PRINT_AND_LOG("Other:" << std::endl); - STORM_PRINT_AND_LOG(" * number of states eliminated: " << maybeStates.getNumberOfSetBits() << std::endl); - STORM_PRINT_AND_LOG(" * maximal depth of SCC decomposition: " << maximalDepth << std::endl); - } - - // Now, we return the value for the only initial state. - return simplify(oneStepProbabilities[initialStateIndex]); + return computeReachabilityProbability(submatrix, oneStepProbabilities, submatrix.transpose(), newInitialStates, phiStates, psiStates, distances); } template @@ -193,7 +174,7 @@ namespace storm { } template - uint_fast64_t SparseSccModelChecker::treatScc(storm::models::Dtmc const& dtmc, FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, std::vector const& distances) { + uint_fast64_t SparseSccModelChecker::treatScc(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, storm::utility::ConstantsComparator const& comparator, boost::optional> const& distances) { uint_fast64_t maximalDepth = level; // If the SCCs are large enough, we try to split them further. @@ -220,7 +201,11 @@ namespace storm { STORM_LOG_DEBUG("Eliminating " << trivialSccs.size() << " trivial SCCs."); if (storm::settings::parametricSettings().isSortTrivialSccsSet()) { - std::sort(trivialSccs.begin(), trivialSccs.end(), [&distances] (std::pair const& state1, std::pair const& state2) -> bool { return distances[state1.first] > distances[state2.first]; } ); + STORM_LOG_THROW(distances, storm::exceptions::IllegalFunctionCallException, "Cannot sort according to distances because none were provided."); + std::vector const& actualDistances = distances.get(); +// std::sort(trivialSccs.begin(), trivialSccs.end(), [&actualDistances] (std::pair const& state1, std::pair const& state2) -> bool { return actualDistances[state1.first] > actualDistances[state2.first]; } ); + + std::sort(trivialSccs.begin(), trivialSccs.end(), [&oneStepProbabilities,&comparator] (std::pair const& state1, std::pair const& state2) -> bool { return comparator.isZero(oneStepProbabilities[state1.first]) && !comparator.isZero(oneStepProbabilities[state2.first]); } ); } for (auto const& stateIndexPair : trivialSccs) { eliminateState(matrix, oneStepProbabilities, stateIndexPair.first, backwardTransitions); @@ -237,7 +222,7 @@ namespace storm { storm::storage::BitVector newSccAsBitVector(forwardTransitions.getRowCount(), newScc.begin(), newScc.end()); // Determine the set of entry states of the SCC. - storm::storage::BitVector entryStates(dtmc.getNumberOfStates()); + storm::storage::BitVector entryStates(forwardTransitions.getRowCount()); for (auto const& state : newScc) { for (auto const& predecessor : backwardTransitions.getRow(state)) { if (predecessor.getValue() != storm::utility::constantZero() && !newSccAsBitVector.get(predecessor.getColumn())) { @@ -247,7 +232,7 @@ namespace storm { } // Recursively descend in SCC-hierarchy. - uint_fast64_t depth = treatScc(dtmc, matrix, oneStepProbabilities, entryStates, newSccAsBitVector, forwardTransitions, backwardTransitions, !storm::settings::parametricSettings().isEliminateEntryStatesLastSet(), level + 1, maximalSccSize, entryStateQueue, distances); + uint_fast64_t depth = treatScc(matrix, oneStepProbabilities, entryStates, newSccAsBitVector, forwardTransitions, backwardTransitions, !storm::settings::parametricSettings().isEliminateEntryStatesLastSet(), level + 1, maximalSccSize, entryStateQueue, comparator, distances); maximalDepth = std::max(maximalDepth, depth); } @@ -256,21 +241,31 @@ namespace storm { STORM_LOG_DEBUG("SCC of size " << scc.getNumberOfSetBits() << " is small enough to be eliminated directly."); storm::storage::BitVector remainingStates = scc & ~entryStates; + std::vector statesToEliminate(remainingStates.begin(), remainingStates.end()); + if (storm::settings::parametricSettings().isSortTrivialSccsSet()) { +// STORM_LOG_THROW(distances, storm::exceptions::IllegalFunctionCallException, "Cannot sort according to distances because none were provided."); +// std::vector const& actualDistances = distances.get(); +// std::sort(statesToEliminate.begin(), statesToEliminate.end(), [&actualDistances] (storm::storage::sparse::state_type const& state1, storm::storage::sparse::state_type const& state2) -> bool { return actualDistances[state1] > actualDistances[state2]; } ); + + std::sort(statesToEliminate.begin(), statesToEliminate.end(), [&oneStepProbabilities,&comparator] (storm::storage::sparse::state_type const& state1, storm::storage::sparse::state_type const& state2) -> bool { return comparator.isZero(oneStepProbabilities[state1]) && !comparator.isZero(oneStepProbabilities[state2]); } ); + + } + // Eliminate the remaining states that do not have a self-loop (in the current, i.e. modified) // transition probability matrix. - for (auto const& state : remainingStates) { - if (!hasSelfLoop(state, matrix)) { + for (auto const& state : statesToEliminate) { +// if (!hasSelfLoop(state, matrix)) { eliminateState(matrix, oneStepProbabilities, state, backwardTransitions); - remainingStates.set(state, false); - } +// remainingStates.set(state, false); +// } } STORM_LOG_DEBUG("Eliminated all states without self-loop."); // Eliminate the remaining states. - for (auto const& state : remainingStates) { - eliminateState(matrix, oneStepProbabilities, state, backwardTransitions); - } +// for (auto const& state : statesToEliminate) { +// eliminateState(matrix, oneStepProbabilities, state, backwardTransitions); +// } STORM_LOG_DEBUG("Eliminated all states with self-loop."); } @@ -335,7 +330,7 @@ namespace storm { std::size_t scaledSuccessors = 0; if (hasSelfLoop) { loopProbability = storm::utility::constantOne() / (storm::utility::constantOne() - loopProbability); - simplify(loopProbability); + storm::utility::simplify(loopProbability); for (auto& entry : matrix.getRow(state)) { // Only scale the non-diagonal entries. if (entry.getColumn() != state) { @@ -409,7 +404,7 @@ namespace storm { auto tmpResult = *first2 * multiplyFactor; multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock; simplifyClock = std::chrono::high_resolution_clock::now(); - *result = simplify(std::move(tmpResult)); + *result = storm::utility::simplify(tmpResult); simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; ++first2; } else if (first1->getColumn() < first2->getColumn()) { @@ -420,14 +415,14 @@ namespace storm { auto tmp1 = multiplyFactor * first2->getValue(); multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock; simplifyClock = std::chrono::high_resolution_clock::now(); - tmp1 = simplify(std::move(tmp1)); + tmp1 = storm::utility::simplify(tmp1); multiplicationClock = std::chrono::high_resolution_clock::now(); simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; additionClock = std::chrono::high_resolution_clock::now(); auto tmp2 = first1->getValue() + tmp1; additionTime += std::chrono::high_resolution_clock::now() - additionClock; simplifyClock = std::chrono::high_resolution_clock::now(); - tmp2 = simplify(std::move(tmp2)); + tmp2 = storm::utility::simplify(tmp2); simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; *result = storm::storage::MatrixEntry::index_type, typename FlexibleSparseMatrix::value_type>(first1->getColumn(), tmp2); ++first1; @@ -440,7 +435,7 @@ namespace storm { auto tmpResult = *first2 * multiplyFactor; multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock; simplifyClock = std::chrono::high_resolution_clock::now(); - *result = simplify(std::move(tmpResult)); + *result = storm::utility::simplify(tmpResult); simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; } } @@ -449,19 +444,19 @@ namespace storm { predecessorForwardTransitions = std::move(newSuccessors); // Add the probabilities to go to a target state in just one step. - multiplicationClock = std::chrono::high_resolution_clock::now(); - auto tmp1 = multiplyFactor * oneStepProbabilities[state]; - multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock; - simplifyClock = std::chrono::high_resolution_clock::now(); - tmp1 = simplify(std::move(tmp1)); - simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; +// multiplicationClock = std::chrono::high_resolution_clock::now(); +// auto tmp1 = multiplyFactor * oneStepProbabilities[state]; +// multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock; +// simplifyClock = std::chrono::high_resolution_clock::now(); +// tmp1 = storm::utility::simplify(tmp1); +// simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; +// auto tmp2 = oneStepProbabilities[predecessor] + tmp1; +// simplifyClock = std::chrono::high_resolution_clock::now(); +// tmp2 = storm::utility::simplify(tmp2); +// simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; additionClock2 = std::chrono::high_resolution_clock::now(); - auto tmp2 = oneStepProbabilities[predecessor] + tmp1; + oneStepProbabilities[predecessor] += storm::utility::simplify(multiplyFactor * oneStepProbabilities[state]); additionTime2 += std::chrono::high_resolution_clock::now() - additionClock2; - simplifyClock = std::chrono::high_resolution_clock::now(); - tmp2 = simplify(std::move(tmp2)); - simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock; - oneStepProbabilities[predecessor] = tmp2; STORM_LOG_DEBUG("Fixed new next-state probabilities of predecessor states."); } @@ -517,6 +512,8 @@ namespace storm { // Clear the eliminated row to reduce memory consumption. currentStateSuccessors.clear(); currentStateSuccessors.shrink_to_fit(); + currentStatePredecessors.clear(); + currentStatePredecessors.shrink_to_fit(); auto eliminationEnd = std::chrono::high_resolution_clock::now(); auto eliminationTime = eliminationEnd - eliminationStart; diff --git a/src/modelchecker/reachability/SparseSccModelChecker.h b/src/modelchecker/reachability/SparseSccModelChecker.h index a12cf9196..039d30cc1 100644 --- a/src/modelchecker/reachability/SparseSccModelChecker.h +++ b/src/modelchecker/reachability/SparseSccModelChecker.h @@ -36,9 +36,11 @@ namespace storm { class SparseSccModelChecker { public: static ValueType computeReachabilityProbability(storm::models::Dtmc const& dtmc, std::shared_ptr> const& filterFormula); - + + static ValueType computeReachabilityProbability(storm::storage::SparseMatrix const& transitionMatrix, std::vector& oneStepProbabilities, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional> const& distances = {}); + private: - static uint_fast64_t treatScc(storm::models::Dtmc const& dtmc, FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, std::vector const& distances); + static uint_fast64_t treatScc(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix const& forwardTransitions, FlexibleSparseMatrix& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector& entryStateQueue, storm::utility::ConstantsComparator const& comparator, boost::optional> const& distances); static FlexibleSparseMatrix getFlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false); static void eliminateState(FlexibleSparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix& backwardTransitions); static bool eliminateStateInPlace(storm::storage::SparseMatrix& matrix, std::vector& oneStepProbabilities, uint_fast64_t state, storm::storage::SparseMatrix& backwardTransitions); diff --git a/src/models/AtomicPropositionsLabeling.h b/src/models/AtomicPropositionsLabeling.h index 1969a525a..dfcc1db04 100644 --- a/src/models/AtomicPropositionsLabeling.h +++ b/src/models/AtomicPropositionsLabeling.h @@ -249,7 +249,7 @@ public: auto apIndexPair = nameToLabelingMap.find(ap); return this->singleLabelings[apIndexPair->second].get(state); } - + /*! * Returns the number of atomic propositions managed by this object (set in the initialization). * diff --git a/src/models/Ctmc.h b/src/models/Ctmc.h index 229ccfe9e..b3c1f1ef1 100644 --- a/src/models/Ctmc.h +++ b/src/models/Ctmc.h @@ -36,8 +36,9 @@ public: * propositions to each state. */ Ctmc(storm::storage::SparseMatrix const& rateMatrix, storm::models::AtomicPropositionsLabeling const& stateLabeling, - boost::optional> const& optionalStateRewardVector, boost::optional> const& optionalTransitionRewardMatrix, - boost::optional>> const& optionalChoiceLabeling) + boost::optional> const& optionalStateRewardVector = boost::optional>(), + boost::optional> const& optionalTransitionRewardMatrix = boost::optional>(), + boost::optional>> const& optionalChoiceLabeling = boost::optional>>()) : AbstractDeterministicModel(rateMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { // Intentionally left empty. } @@ -51,8 +52,9 @@ public: * propositions to each state. */ Ctmc(storm::storage::SparseMatrix&& rateMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling, - boost::optional>&& optionalStateRewardVector, boost::optional>&& optionalTransitionRewardMatrix, - boost::optional>>&& optionalChoiceLabeling) + boost::optional>&& optionalStateRewardVector = boost::optional>(), + boost::optional>&& optionalTransitionRewardMatrix = boost::optional>(), + boost::optional>>&& optionalChoiceLabeling = boost::optional>>()) // The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class : AbstractDeterministicModel(std::move(rateMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h index dcb804fe8..d3d4a1b8f 100644 --- a/src/models/Dtmc.h +++ b/src/models/Dtmc.h @@ -20,6 +20,7 @@ #include "src/settings/SettingsManager.h" #include "src/utility/vector.h" #include "src/utility/matrix.h" +#include "src/utility/ConstantsComparator.h" namespace storm { @@ -73,8 +74,9 @@ public: * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state. */ Dtmc(storm::storage::SparseMatrix&& probabilityMatrix, storm::models::AtomicPropositionsLabeling&& stateLabeling, - boost::optional>&& optionalStateRewardVector, boost::optional>&& optionalTransitionRewardMatrix, - boost::optional>>&& optionalChoiceLabeling) + boost::optional>&& optionalStateRewardVector = boost::optional>(), + boost::optional>&& optionalTransitionRewardMatrix = boost::optional>(), + boost::optional>>&& optionalChoiceLabeling = boost::optional>>()) // The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class : AbstractDeterministicModel(std::move(probabilityMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { @@ -328,24 +330,23 @@ private: * Checks probability matrix if all rows sum up to one. */ bool checkValidityOfProbabilityMatrix() { - // Get the settings object to customize linear solving. - - if (this->getTransitionMatrix().getRowCount() != this->getTransitionMatrix().getColumnCount()) { // not square LOG4CPLUS_ERROR(logger, "Probability matrix is not square."); return false; } + + storm::utility::ConstantsComparator comparator; for (uint_fast64_t row = 0; row < this->getTransitionMatrix().getRowCount(); ++row) { - T sum = this->getTransitionMatrix().getRowSum(row); - - if (sum == T(0)) { - - LOG4CPLUS_ERROR(logger, "Row " << row << " is a deadlock (sum == " << sum << ")."); - return false; - } - if (!storm::utility::isOne(sum)) { - LOG4CPLUS_ERROR(logger, "Row " << row << " has sum " << sum << "."); + T sum = this->getTransitionMatrix().getRowSum(row); + + // If the sum is not a constant, for example for parametric models, we cannot check whether the sum is one + // or not. + if (!comparator.isConstant(sum)) { + continue; + } + + if (!comparator.isOne(sum)) { return false; } } diff --git a/src/parser/AutoParser.cpp b/src/parser/AutoParser.cpp index 697740a7f..fbcc0b393 100644 --- a/src/parser/AutoParser.cpp +++ b/src/parser/AutoParser.cpp @@ -85,7 +85,7 @@ namespace storm { // Find and read in the hint. std::string formatString = "%" + std::to_string(STORM_PARSER_AUTOPARSER_HINT_LENGTH) + "s"; - char hint[5]; + char hint[STORM_PARSER_AUTOPARSER_HINT_LENGTH + 1]; #ifdef WINDOWS sscanf_s(filehintBuffer, formatString.c_str(), hint, STORM_PARSER_AUTOPARSER_HINT_LENGTH + 1); #else diff --git a/src/parser/CslParser.cpp b/src/parser/CslParser.cpp index 4aa8f7d91..b16e441e8 100644 --- a/src/parser/CslParser.cpp +++ b/src/parser/CslParser.cpp @@ -80,8 +80,7 @@ struct CslParser::CslGrammar : qi::grammar> stateFormula >> qi::lit(")")) | (qi::lit("[") >> stateFormula >> qi::lit("]")); atomicStateFormula.name("atomic state formula"); - atomicProposition = (freeIdentifierName)[qi::_val = - MAKE(csl::Ap, qi::_1)]; + atomicProposition = (freeIdentifierName)[qi::_val = MAKE(csl::Ap, qi::_1)] | (qi::lit("\"") > (freeIdentifierName)[qi::_val = MAKE(csl::Ap, qi::_1)] > qi::lit("\"")); atomicProposition.name("atomic proposition"); probabilisticBoundOperator = ( (qi::lit("P") >> comparisonType > qi::double_ > pathFormula )[qi::_val = diff --git a/src/parser/LtlParser.cpp b/src/parser/LtlParser.cpp index 86fff3b03..79059e994 100644 --- a/src/parser/LtlParser.cpp +++ b/src/parser/LtlParser.cpp @@ -88,8 +88,7 @@ struct LtlParser::LtlGrammar : qi::grammar> formula >> qi::lit(")")| qi::lit("[") >> formula >> qi::lit("]"); atomicLtlFormula.name("Atomic LTL formula"); - atomicProposition = (freeIdentifierName)[qi::_val = - MAKE(ltl::Ap, qi::_1)]; + atomicProposition = (freeIdentifierName)[qi::_val = MAKE(ltl::Ap, qi::_1)] | (qi::lit("\"") > (freeIdentifierName)[qi::_val = MAKE(ltl::Ap, qi::_1)] > qi::lit("\"")); atomicProposition.name("Atomic Proposition"); //This block defines rules for parsing probabilistic path formulas diff --git a/src/parser/PrctlParser.cpp b/src/parser/PrctlParser.cpp index d2e6f06fa..d12fba12b 100644 --- a/src/parser/PrctlParser.cpp +++ b/src/parser/PrctlParser.cpp @@ -76,8 +76,7 @@ struct PrctlParser::PrctlGrammar : qi::grammar> stateFormula >> qi::lit(")") | qi::lit("[") >> stateFormula >> qi::lit("]"); atomicStateFormula.name("atomic state formula"); - atomicProposition = (freeIdentifierName)[qi::_val = - MAKE(prctl::Ap, qi::_1)]; + atomicProposition = (freeIdentifierName)[qi::_val = MAKE(prctl::Ap, qi::_1)] | (qi::lit("\"") > (freeIdentifierName)[qi::_val = MAKE(prctl::Ap, qi::_1)] > qi::lit("\"")); atomicProposition.name("atomic proposition"); probabilisticBoundOperator = ((qi::lit("P") >> comparisonType > qi::double_ > pathFormula)[qi::_val = MAKE(prctl::ProbabilisticBoundOperator, qi::_1, qi::_2, qi::_3)]); diff --git a/src/properties/prctl/PrctlFilter.h b/src/properties/prctl/PrctlFilter.h index dc7c2c382..e72d700c8 100644 --- a/src/properties/prctl/PrctlFilter.h +++ b/src/properties/prctl/PrctlFilter.h @@ -194,7 +194,6 @@ namespace storm { } else { if(this->actions.empty()) { - // There are no filter actions but only the raw state formula. So just print that. return child->toString(); } diff --git a/src/settings/modules/GeneralSettings.cpp b/src/settings/modules/GeneralSettings.cpp index 0e1b8775a..d04ec0d25 100644 --- a/src/settings/modules/GeneralSettings.cpp +++ b/src/settings/modules/GeneralSettings.cpp @@ -43,6 +43,8 @@ namespace storm { const std::string GeneralSettings::constantsOptionShortName = "const"; const std::string GeneralSettings::statisticsOptionName = "statistics"; const std::string GeneralSettings::statisticsOptionShortName = "stats"; + const std::string GeneralSettings::bisimulationOptionName = "bisimulation"; + const std::string GeneralSettings::bisimulationOptionShortName = "bisim"; GeneralSettings::GeneralSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) { this->addOption(storm::settings::OptionBuilder(moduleName, helpOptionName, false, "Shows all available options, arguments and descriptions.").setShortName(helpOptionShortName) @@ -74,7 +76,7 @@ namespace storm { .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the LTL formulas.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, counterexampleOptionName, false, "Generates a counterexample for the given PRCTL formulas if not satisfied by the model") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file to which the counterexample is to be written.").setDefaultValueString("-").setIsOptional(true).build()).setShortName(counterexampleOptionShortName).build()); - + this->addOption(storm::settings::OptionBuilder(moduleName, bisimulationOptionName, false, "Sets whether to perform bisimulation minimization.").setShortName(bisimulationOptionShortName).build()); this->addOption(storm::settings::OptionBuilder(moduleName, transitionRewardsOptionName, "", "If given, the transition rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the transition rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, stateRewardsOptionName, false, "If given, the state rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").") @@ -93,7 +95,6 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, constantsOptionName, false, "Specifies the constant replacements to use in symbolic models. Note that Note that this requires the model to be given as an symbolic model (i.e., via --" + symbolicOptionName + ").").setShortName(constantsOptionShortName) .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of constants and their value, e.g. a=1,b=2,c=3.").setDefaultValueString("").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, statisticsOptionName, false, "Sets whether to display statistics if available.").setShortName(statisticsOptionShortName).build()); - } bool GeneralSettings::isHelpSet() const { @@ -284,6 +285,10 @@ namespace storm { return true; } + bool GeneralSettings::isBisimulationSet() const { + return this->getOption(bisimulationOptionName).getHasOptionBeenSet(); + } + } // namespace modules } // namespace settings } // namespace storm \ No newline at end of file diff --git a/src/settings/modules/GeneralSettings.h b/src/settings/modules/GeneralSettings.h index 912053083..c65ab57a2 100644 --- a/src/settings/modules/GeneralSettings.h +++ b/src/settings/modules/GeneralSettings.h @@ -322,6 +322,13 @@ namespace storm { */ bool isShowStatisticsSet() const; + /*! + * Retrieves whether the option to perform bisimulation minimization is set. + * + * @return True iff the option was set. + */ + bool isBisimulationSet() const; + bool check() const override; // The name of the module. @@ -363,6 +370,8 @@ namespace storm { static const std::string constantsOptionShortName; static const std::string statisticsOptionName; static const std::string statisticsOptionShortName; + static const std::string bisimulationOptionName; + static const std::string bisimulationOptionShortName; }; } // namespace modules diff --git a/src/storage/Decomposition.cpp b/src/storage/Decomposition.cpp index 06656b483..39c49bf4f 100644 --- a/src/storage/Decomposition.cpp +++ b/src/storage/Decomposition.cpp @@ -89,7 +89,10 @@ namespace storm { out << "]"; return out; } - + + template class Decomposition; + template std::ostream& operator<<(std::ostream& out, Decomposition const& decomposition); + template class Decomposition; template std::ostream& operator<<(std::ostream& out, Decomposition const& decomposition); diff --git a/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp new file mode 100644 index 000000000..1d455be2d --- /dev/null +++ b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp @@ -0,0 +1,943 @@ +#include "src/storage/DeterministicModelStrongBisimulationDecomposition.h" + +#include +#include +#include +#include + +#include "src/utility/graph.h" +#include "src/exceptions/IllegalFunctionCallException.h" + +namespace storm { + namespace storage { + + template + std::size_t DeterministicModelStrongBisimulationDecomposition::Block::blockId = 0; + + template + DeterministicModelStrongBisimulationDecomposition::Block::Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::shared_ptr const& label) : next(next), prev(prev), begin(begin), end(end), markedAsSplitter(false), markedAsPredecessorBlock(false), markedPosition(begin), absorbing(false), id(blockId++), label(label) { + if (next != nullptr) { + next->prev = this; + } + if (prev != nullptr) { + prev->next = this; + } + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::print(Partition const& partition) const { + std::cout << "block " << this->getId() << " with marked position " << this->getMarkedPosition() << " and original begin " << this->getOriginalBegin() << std::endl; + std::cout << "begin: " << this->begin << " and end: " << this->end << " (number of states: " << this->getNumberOfStates() << ")" << std::endl; + std::cout << "states:" << std::endl; + for (storm::storage::sparse::state_type index = this->begin; index < this->end; ++index) { + std::cout << partition.states[index] << " "; + } + std::cout << std::endl << "original: " << std::endl; + for (storm::storage::sparse::state_type index = this->getOriginalBegin(); index < this->end; ++index) { + std::cout << partition.states[index] << " "; + } + std::cout << std::endl << "values:" << std::endl; + for (storm::storage::sparse::state_type index = this->begin; index < this->end; ++index) { + std::cout << std::setprecision(3) << partition.values[index] << " "; + } + std::cout << std::endl; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::setBegin(storm::storage::sparse::state_type begin) { + this->begin = begin; + this->markedPosition = begin; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::setEnd(storm::storage::sparse::state_type end) { + this->end = end; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::incrementBegin() { + ++this->begin; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::decrementEnd() { + ++this->begin; + } + + template + storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition::Block::getBegin() const { + return this->begin; + } + + template + storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition::Block::getOriginalBegin() const { + if (this->hasPreviousBlock()) { + return this->getPreviousBlock().getEnd(); + } else { + return 0; + } + } + + template + storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition::Block::getEnd() const { + return this->end; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::setIterator(const_iterator it) { + this->selfIt = it; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block::const_iterator DeterministicModelStrongBisimulationDecomposition::Block::getIterator() const { + return this->selfIt; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block::const_iterator DeterministicModelStrongBisimulationDecomposition::Block::getNextIterator() const { + return this->getNextBlock().getIterator(); + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block::const_iterator DeterministicModelStrongBisimulationDecomposition::Block::getPreviousIterator() const { + return this->getPreviousBlock().getIterator(); + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Block::getNextBlock() { + return *this->next; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block const& DeterministicModelStrongBisimulationDecomposition::Block::getNextBlock() const { + return *this->next; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::hasNextBlock() const { + return this->next != nullptr; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Block::getPreviousBlock() { + return *this->prev; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block* DeterministicModelStrongBisimulationDecomposition::Block::getPreviousBlockPointer() { + return this->prev; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block* DeterministicModelStrongBisimulationDecomposition::Block::getNextBlockPointer() { + return this->next; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block const& DeterministicModelStrongBisimulationDecomposition::Block::getPreviousBlock() const { + return *this->prev; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::hasPreviousBlock() const { + return this->prev != nullptr; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::check() const { + if (this->begin >= this->end) { + assert(false); + } + if (this->prev != nullptr) { + assert (this->prev->next == this); + } + if (this->next != nullptr) { + assert (this->next->prev == this); + } + return true; + } + + template + std::size_t DeterministicModelStrongBisimulationDecomposition::Block::getNumberOfStates() const { + // We need to take the original begin here, because the begin is temporarily moved. + return (this->end - this->getOriginalBegin()); + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::isMarkedAsSplitter() const { + return this->markedAsSplitter; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::markAsSplitter() { + this->markedAsSplitter = true; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::unmarkAsSplitter() { + this->markedAsSplitter = false; + } + + template + std::size_t DeterministicModelStrongBisimulationDecomposition::Block::getId() const { + return this->id; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::setAbsorbing(bool absorbing) { + this->absorbing = absorbing; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::isAbsorbing() const { + return this->absorbing; + } + + template + storm::storage::sparse::state_type DeterministicModelStrongBisimulationDecomposition::Block::getMarkedPosition() const { + return this->markedPosition; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::setMarkedPosition(storm::storage::sparse::state_type position) { + this->markedPosition = position; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::resetMarkedPosition() { + this->markedPosition = this->begin; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::incrementMarkedPosition() { + ++this->markedPosition; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::markAsPredecessorBlock() { + this->markedAsPredecessorBlock = true; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Block::unmarkAsPredecessorBlock() { + this->markedAsPredecessorBlock = false; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::isMarkedAsPredecessor() const { + return markedAsPredecessorBlock; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Block::hasLabel() const { + return this->label != nullptr; + } + + template + std::string const& DeterministicModelStrongBisimulationDecomposition::Block::getLabel() const { + STORM_LOG_THROW(this->label != nullptr, storm::exceptions::IllegalFunctionCallException, "Unable to retrieve label of block that has none."); + return *this->label; + } + + template + std::shared_ptr const& DeterministicModelStrongBisimulationDecomposition::Block::getLabelPtr() const { + return this->label; + } + + template + DeterministicModelStrongBisimulationDecomposition::Partition::Partition(std::size_t numberOfStates) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates), values(numberOfStates) { + // Create the block and give it an iterator to itself. + typename std::list::iterator it = blocks.emplace(this->blocks.end(), 0, numberOfStates, nullptr, nullptr); + it->setIterator(it); + + // Set up the different parts of the internal structure. + for (storm::storage::sparse::state_type state = 0; state < numberOfStates; ++state) { + states[state] = state; + positions[state] = state; + stateToBlockMapping[state] = &blocks.back(); + } + } + + template + DeterministicModelStrongBisimulationDecomposition::Partition::Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, std::string const& otherLabel, std::string const& prob1Label) : stateToBlockMapping(numberOfStates), states(numberOfStates), positions(numberOfStates), values(numberOfStates) { + typename std::list::iterator firstIt = blocks.emplace(this->blocks.end(), 0, prob0States.getNumberOfSetBits(), nullptr, nullptr); + Block& firstBlock = *firstIt; + firstBlock.setIterator(firstIt); + + storm::storage::sparse::state_type position = 0; + for (auto state : prob0States) { + states[position] = state; + positions[state] = position; + stateToBlockMapping[state] = &firstBlock; + ++position; + } + firstBlock.setAbsorbing(true); + + typename std::list::iterator secondIt = blocks.emplace(this->blocks.end(), position, position + prob1States.getNumberOfSetBits(), &firstBlock, nullptr, std::shared_ptr(new std::string(prob1Label))); + Block& secondBlock = *secondIt; + secondBlock.setIterator(secondIt); + + for (auto state : prob1States) { + states[position] = state; + positions[state] = position; + stateToBlockMapping[state] = &secondBlock; + ++position; + } + secondBlock.setAbsorbing(true); + + typename std::list::iterator thirdIt = blocks.emplace(this->blocks.end(), position, numberOfStates, &secondBlock, nullptr, otherLabel == "true" ? std::shared_ptr(nullptr) : std::shared_ptr(new std::string(otherLabel))); + Block& thirdBlock = *thirdIt; + thirdBlock.setIterator(thirdIt); + + storm::storage::BitVector otherStates = ~(prob0States | prob1States); + for (auto state : otherStates) { + states[position] = state; + positions[state] = position; + stateToBlockMapping[state] = &thirdBlock; + ++position; + } + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { + std::swap(this->states[this->positions[state1]], this->states[this->positions[state2]]); + std::swap(this->values[this->positions[state1]], this->values[this->positions[state2]]); + std::swap(this->positions[state1], this->positions[state2]); + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2) { + storm::storage::sparse::state_type state1 = this->states[position1]; + storm::storage::sparse::state_type state2 = this->states[position2]; + + std::swap(this->states[position1], this->states[position2]); + std::swap(this->values[position1], this->values[position2]); + + this->positions[state1] = position2; + this->positions[state2] = position1; + } + + template + storm::storage::sparse::state_type const& DeterministicModelStrongBisimulationDecomposition::Partition::getPosition(storm::storage::sparse::state_type state) const { + return this->positions[state]; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position) { + this->positions[state] = position; + } + + template + storm::storage::sparse::state_type const& DeterministicModelStrongBisimulationDecomposition::Partition::getState(storm::storage::sparse::state_type position) const { + return this->states[position]; + } + + template + ValueType const& DeterministicModelStrongBisimulationDecomposition::Partition::getValue(storm::storage::sparse::state_type state) const { + return this->values[this->positions[state]]; + } + + template + ValueType const& DeterministicModelStrongBisimulationDecomposition::Partition::getValueAtPosition(storm::storage::sparse::state_type position) const { + return this->values[position]; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::setValue(storm::storage::sparse::state_type state, ValueType value) { + this->values[this->positions[state]] = value; + } + + template + std::vector& DeterministicModelStrongBisimulationDecomposition::Partition::getValues() { + return this->values; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::increaseValue(storm::storage::sparse::state_type state, ValueType value) { + this->values[this->positions[state]] += value; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::updateBlockMapping(Block& block, std::vector::iterator first, std::vector::iterator last) { + for (; first != last; ++first) { + this->stateToBlockMapping[*first] = █ + } + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Partition::getFirstBlock() { + return *this->blocks.begin(); + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Partition::getBlock(storm::storage::sparse::state_type state) { + return *this->stateToBlockMapping[state]; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block const& DeterministicModelStrongBisimulationDecomposition::Partition::getBlock(storm::storage::sparse::state_type state) const { + return *this->stateToBlockMapping[state]; + } + + template + std::vector::iterator DeterministicModelStrongBisimulationDecomposition::Partition::getBeginOfStates(Block const& block) { + return this->states.begin() + block.getBegin(); + } + + template + std::vector::iterator DeterministicModelStrongBisimulationDecomposition::Partition::getEndOfStates(Block const& block) { + return this->states.begin() + block.getEnd(); + } + + template + std::vector::const_iterator DeterministicModelStrongBisimulationDecomposition::Partition::getBeginOfStates(Block const& block) const { + return this->states.begin() + block.getBegin(); + } + + template + std::vector::const_iterator DeterministicModelStrongBisimulationDecomposition::Partition::getEndOfStates(Block const& block) const { + return this->states.begin() + block.getEnd(); + } + + template + typename std::vector::iterator DeterministicModelStrongBisimulationDecomposition::Partition::getBeginOfValues(Block const& block) { + return this->values.begin() + block.getBegin(); + } + + template + typename std::vector::iterator DeterministicModelStrongBisimulationDecomposition::Partition::getEndOfValues(Block const& block) { + return this->values.begin() + block.getEnd(); + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Partition::splitBlock(Block& block, storm::storage::sparse::state_type position) { + // FIXME: this could be improved by splitting off the smaller of the two resulting blocks. + + // In case one of the resulting blocks would be empty, we simply return the current block and do not create + // a new one. + if (position == block.getBegin() || position == block.getEnd()) { + return block; + } + + // Actually create the new block and insert it at the correct position. + typename std::list::iterator selfIt = this->blocks.emplace(block.getIterator(), block.getBegin(), position, block.getPreviousBlockPointer(), &block, block.getLabelPtr()); + selfIt->setIterator(selfIt); + Block& newBlock = *selfIt; + + // Make the current block end at the given position. + block.setBegin(position); + + // Update the mapping of the states in the newly created block. + for (auto it = this->getBeginOfStates(newBlock), ite = this->getEndOfStates(newBlock); it != ite; ++it) { + stateToBlockMapping[*it] = &newBlock; + } + + return newBlock; + } + + template + typename DeterministicModelStrongBisimulationDecomposition::Block& DeterministicModelStrongBisimulationDecomposition::Partition::insertBlock(Block& block) { + // Find the beginning of the new block. + storm::storage::sparse::state_type begin; + if (block.hasPreviousBlock()) { + begin = block.getPreviousBlock().getEnd(); + } else { + begin = 0; + } + + // Actually insert the new block. + typename std::list::iterator it = this->blocks.emplace(block.getIterator(), begin, block.getBegin(), block.getPreviousBlockPointer(), &block); + Block& newBlock = *it; + newBlock.setIterator(it); + + // Update the mapping of the states in the newly created block. + for (auto it = this->getBeginOfStates(newBlock), ite = this->getEndOfStates(newBlock); it != ite; ++it) { + stateToBlockMapping[*it] = &newBlock; + } + + return *it; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::splitLabel(storm::storage::BitVector const& statesWithLabel) { + for (auto blockIterator = this->blocks.begin(), ite = this->blocks.end(); blockIterator != ite; ) { // The update of the loop was intentionally moved to the bottom of the loop. + Block& block = *blockIterator; + + // Sort the range of the block such that all states that have the label are moved to the front. + std::sort(this->getBeginOfStates(block), this->getEndOfStates(block), [&statesWithLabel] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return statesWithLabel.get(a) && !statesWithLabel.get(b); } ); + + // Update the positions vector. + storm::storage::sparse::state_type position = block.getBegin(); + for (auto stateIt = this->getBeginOfStates(block), stateIte = this->getEndOfStates(block); stateIt != stateIte; ++stateIt, ++position) { + this->positions[*stateIt] = position; + } + + // Now we can find the first position in the block that does not have the label and create new blocks. + std::vector::iterator it = std::find_if(this->getBeginOfStates(block), this->getEndOfStates(block), [&] (storm::storage::sparse::state_type const& a) { return !statesWithLabel.get(a); }); + + // If not all the states agreed on the validity of the label, we need to split the block. + if (it != this->getBeginOfStates(block) && it != this->getEndOfStates(block)) { + auto cutPoint = std::distance(this->states.begin(), it); + this->splitBlock(block, cutPoint); + } else { + // Otherwise, we simply proceed to the next block. + ++blockIterator; + } + } + } + + template + std::list::Block> const& DeterministicModelStrongBisimulationDecomposition::Partition::getBlocks() const { + return this->blocks; + } + + template + std::vector& DeterministicModelStrongBisimulationDecomposition::Partition::getStates() { + return this->states; + } + + template + std::list::Block>& DeterministicModelStrongBisimulationDecomposition::Partition::getBlocks() { + return this->blocks; + } + + template + bool DeterministicModelStrongBisimulationDecomposition::Partition::check() const { + for (uint_fast64_t state = 0; state < this->positions.size(); ++state) { + if (this->states[this->positions[state]] != state) { + assert(false); + } + } + for (auto const& block : this->blocks) { + assert(block.check()); + for (auto stateIt = this->getBeginOfStates(block), stateIte = this->getEndOfStates(block); stateIt != stateIte; ++stateIt) { + if (this->stateToBlockMapping[*stateIt] != &block) { + assert(false); + } + } + } + return true; + } + + template + void DeterministicModelStrongBisimulationDecomposition::Partition::print() const { + for (auto const& block : this->blocks) { + block.print(*this); + } + std::cout << "states in partition" << std::endl; + for (auto const& state : states) { + std::cout << state << " "; + } + std::cout << std::endl << "positions: " << std::endl; + for (auto const& index : positions) { + std::cout << index << " "; + } + std::cout << std::endl << "state to block mapping: " << std::endl; + for (auto const& block : stateToBlockMapping) { + std::cout << block << " "; + } + std::cout << std::endl; + std::cout << "size: " << this->blocks.size() << std::endl; + assert(this->check()); + } + + template + std::size_t DeterministicModelStrongBisimulationDecomposition::Partition::size() const { + return this->blocks.size(); + } + + template + DeterministicModelStrongBisimulationDecomposition::DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc const& model, bool buildQuotient) { + STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures."); + Partition initialPartition = getLabelBasedInitialPartition(model); + partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient); + } + + template + DeterministicModelStrongBisimulationDecomposition::DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc const& model, bool buildQuotient) { + STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures."); + Partition initialPartition = getLabelBasedInitialPartition(model); + partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient); + } + + template + DeterministicModelStrongBisimulationDecomposition::DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient) { + STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures."); + storm::storage::SparseMatrix backwardTransitions = model.getBackwardTransitions(); + Partition initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, phiLabel, psiLabel, bounded); + partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient); + } + + template + DeterministicModelStrongBisimulationDecomposition::DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient) { + STORM_LOG_THROW(!model.hasStateRewards() && !model.hasTransitionRewards(), storm::exceptions::IllegalFunctionCallException, "Bisimulation is currently only supported for models without reward structures."); + storm::storage::SparseMatrix backwardTransitions = model.getBackwardTransitions(); + Partition initialPartition = getMeasureDrivenInitialPartition(model, backwardTransitions, phiLabel, psiLabel, bounded); + partitionRefinement(model, model.getBackwardTransitions(), initialPartition, buildQuotient); + } + + template + std::shared_ptr> DeterministicModelStrongBisimulationDecomposition::getQuotient() const { + STORM_LOG_THROW(this->quotient != nullptr, storm::exceptions::IllegalFunctionCallException, "Unable to retrieve quotient model from bisimulation decomposition, because it was not built."); + return this->quotient; + } + + template + template + void DeterministicModelStrongBisimulationDecomposition::buildQuotient(ModelType const& model, Partition const& partition) { + // In order to create the quotient model, we need to construct + // (a) the new transition matrix, + // (b) the new labeling, + // (c) the new reward structures. + + // Prepare a matrix builder for (a). + storm::storage::SparseMatrixBuilder builder(this->size(), this->size()); + + // Prepare the new state labeling for (b). + storm::models::AtomicPropositionsLabeling newLabeling(this->size(), model.getStateLabeling().getNumberOfAtomicPropositions()); + std::set atomicPropositionsSet = model.getStateLabeling().getAtomicPropositions(); + std::vector atomicPropositions = std::vector(atomicPropositionsSet.begin(), atomicPropositionsSet.end()); + for (auto const& ap : atomicPropositions) { + newLabeling.addAtomicProposition(ap); + } + + // Now build (a) and (b) by traversing all blocks. + for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) { + auto const& block = this->blocks[blockIndex]; + + // Pick one representative state. It doesn't matter which state it is, because they all behave equally. + storm::storage::sparse::state_type representativeState = *block.begin(); + Block const& oldBlock = partition.getBlock(representativeState); + + // If the block is absorbing, we simply add a self-loop. + if (oldBlock.isAbsorbing()) { + builder.addNextValue(blockIndex, blockIndex, storm::utility::constantOne()); + + if (oldBlock.hasLabel()) { + newLabeling.addAtomicPropositionToState(oldBlock.getLabel(), blockIndex); + } + } else { + // Compute the outgoing transitions of the block. + std::map blockProbability; + for (auto const& entry : model.getTransitionMatrix().getRow(representativeState)) { + storm::storage::sparse::state_type targetBlock = partition.getBlock(entry.getColumn()).getId(); + auto probIterator = blockProbability.find(targetBlock); + if (probIterator != blockProbability.end()) { + probIterator->second += entry.getValue(); + } else { + blockProbability[targetBlock] = entry.getValue(); + } + } + +#ifdef DEBUG + // FIXME: take this out. + // As a debug check, we walk through all states and check if their behaviour is the same modulo the + // partition. + for (auto const& state : block) { + std::map stateProbability; + for (auto const& entry : model.getTransitionMatrix().getRow(state)) { + storm::storage::sparse::state_type targetBlock = partition.getBlock(entry.getColumn()).getId(); + auto probIterator = stateProbability.find(targetBlock); + if (probIterator != stateProbability.end()) { + probIterator->second += entry.getValue(); + } else { + stateProbability[targetBlock] = entry.getValue(); + } + } + if (stateProbability != blockProbability) { + std::cout << "state: " << std::endl; + for (auto const& entry : stateProbability) { + std::cout << entry.first << " -> " << entry.second << std::endl; + } + std::cout << "block: " << std::endl; + for (auto const& entry : blockProbability) { + std::cout << entry.first << " -> " << entry.second << std::endl; + } + } + STORM_LOG_ASSERT(stateProbability == blockProbability, "Quotient distributions did not match."); + } +#endif + + // Now add them to the actual matrix. + for (auto const& probabilityEntry : blockProbability) { + builder.addNextValue(blockIndex, probabilityEntry.first, probabilityEntry.second); + } + + // If the block has a special label, we only add that to the block. + if (oldBlock.hasLabel()) { + newLabeling.addAtomicPropositionToState(oldBlock.getLabel(), blockIndex); + } else { + // Otherwise add all atomic propositions to the equivalence class that the representative state + // satisfies. + for (auto const& ap : atomicPropositions) { + if (model.getStateLabeling().getStateHasAtomicProposition(ap, representativeState)) { + newLabeling.addAtomicPropositionToState(ap, blockIndex); + } + } + } + } + } + + // Now check which of the blocks of the partition contain at least one initial state. + for (auto initialState : model.getInitialStates()) { + Block const& initialBlock = partition.getBlock(initialState); + newLabeling.addAtomicPropositionToState("init", initialBlock.getId()); + } + + // FIXME: + // If reward structures are allowed, the quotient structures need to be built here. + + // Finally construct the quotient model. + this->quotient = std::shared_ptr>(new ModelType(builder.build(), std::move(newLabeling))); + } + + template + template + void DeterministicModelStrongBisimulationDecomposition::partitionRefinement(ModelType const& model, storm::storage::SparseMatrix const& backwardTransitions, Partition& partition, bool buildQuotient) { + std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now(); + + // Initially, all blocks are potential splitter, so we insert them in the splitterQueue. + std::deque splitterQueue; + std::for_each(partition.getBlocks().begin(), partition.getBlocks().end(), [&] (Block& a) { splitterQueue.push_back(&a); }); + + // Then perform the actual splitting until there are no more splitters. + while (!splitterQueue.empty()) { + refinePartition(backwardTransitions, *splitterQueue.front(), partition, splitterQueue); + splitterQueue.pop_front(); + } + + // Now move the states from the internal partition into their final place in the decomposition. We do so in + // a way that maintains the block IDs as indices. + this->blocks.resize(partition.size()); + for (auto const& block : partition.getBlocks()) { + // We need to sort the states to allow for rapid construction of the blocks. + std::sort(partition.getBeginOfStates(block), partition.getEndOfStates(block)); + this->blocks[block.getId()] = block_type(partition.getBeginOfStates(block), partition.getEndOfStates(block), true); + } + + // If we are required to build the quotient model, do so now. + if (buildQuotient) { + this->buildQuotient(model, partition); + } + + std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart; + + if (storm::settings::generalSettings().isShowStatisticsSet()) { + std::cout << "Computed bisimulation quotient in " << std::chrono::duration_cast(totalTime).count() << "ms." << std::endl; + } + } + + template + void DeterministicModelStrongBisimulationDecomposition::refineBlockProbabilities(Block& block, Partition& partition, std::deque& splitterQueue) { + // First, we simplify all the values. For most types this does not change anything, but for rational + // functions, this achieves a canonicity that is used for sorting the rational functions later. + for (auto probabilityIt = partition.getBeginOfValues(block), probabilityIte = partition.getEndOfValues(block); probabilityIt != probabilityIte; ++probabilityIt) { + storm::utility::simplify(*probabilityIt); + } + + // Sort the states in the block based on their probabilities. + std::sort(partition.getBeginOfStates(block), partition.getEndOfStates(block), [&partition] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return partition.getValue(a) < partition.getValue(b); } ); + + // FIXME: This can probably be done more efficiently. + std::sort(partition.getBeginOfValues(block), partition.getEndOfValues(block)); + + // Update the positions vector. + storm::storage::sparse::state_type position = block.getBegin(); + for (auto stateIt = partition.getBeginOfStates(block), stateIte = partition.getEndOfStates(block); stateIt != stateIte; ++stateIt, ++position) { + partition.setPosition(*stateIt, position); + } + + // Finally, we need to scan the ranges of states that agree on the probability. + storm::storage::sparse::state_type beginIndex = block.getBegin(); + storm::storage::sparse::state_type currentIndex = beginIndex; + storm::storage::sparse::state_type endIndex = block.getEnd(); + + // Now we can check whether the block needs to be split, which is the case iff the probabilities for the + // first and the last state are different. + while (!comparator.isEqual(partition.getValueAtPosition(beginIndex), partition.getValueAtPosition(endIndex - 1))) { + // Now we scan for the first state in the block that disagrees on the probability value. + // Note that we do not have to check currentIndex for staying within bounds, because we know the matching + // state is within bounds. + ValueType const& currentValue = partition.getValueAtPosition(beginIndex); + + typename std::vector::const_iterator valueIterator = partition.getValues().begin() + beginIndex + 1; + ++currentIndex; + while (currentIndex < endIndex && comparator.isEqual(*valueIterator, currentValue)) { + ++valueIterator; + ++currentIndex; + } + + // Now we split the block and mark it as a potential splitter. + Block& newBlock = partition.splitBlock(block, currentIndex); + if (!newBlock.isMarkedAsSplitter()) { + splitterQueue.push_back(&newBlock); + newBlock.markAsSplitter(); + } + beginIndex = currentIndex; + } + + for (auto index = currentIndex; index < endIndex - 1; ++index) { + STORM_LOG_ASSERT(comparator.isEqual(partition.getValueAtPosition(index + 1), partition.getValueAtPosition(index + 1)), "Expected equal probabilities after sorting block, but got '" << partition.getValueAtPosition(index) << "' and '" << partition.getValueAtPosition(index + 1) << "'."); + } + } + + template + void DeterministicModelStrongBisimulationDecomposition::refinePartition(storm::storage::SparseMatrix const& backwardTransitions, Block& splitter, Partition& partition, std::deque& splitterQueue) { + std::list predecessorBlocks; + + // Iterate over all states of the splitter and check its predecessors. + storm::storage::sparse::state_type currentPosition = splitter.getBegin(); + for (auto stateIterator = partition.getBeginOfStates(splitter), stateIte = partition.getEndOfStates(splitter); stateIterator != stateIte; ++stateIterator, ++currentPosition) { + storm::storage::sparse::state_type currentState = *stateIterator; + + uint_fast64_t elementsToSkip = 0; + for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) { + storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn(); + + Block& predecessorBlock = partition.getBlock(predecessor); + + // If the predecessor block has just one state, there is no point in splitting it. + if (predecessorBlock.getNumberOfStates() <= 1 || predecessorBlock.isAbsorbing()) { + continue; + } + + storm::storage::sparse::state_type predecessorPosition = partition.getPosition(predecessor); + + // If we have not seen this predecessor before, we move it to a part near the beginning of the block. + if (predecessorPosition >= predecessorBlock.getBegin()) { + if (&predecessorBlock == &splitter) { + // If the predecessor we just found was already processed (in terms of visiting its predecessors), + // we swap it with the state that is currently at the beginning of the block and move the start + // of the block one step further. + if (predecessorPosition <= currentPosition + elementsToSkip) { + partition.swapStates(predecessor, partition.getState(predecessorBlock.getBegin())); + predecessorBlock.incrementBegin(); + } else { + // Otherwise, we need to move the predecessor, but we need to make sure that we explore its + // predecessors later. + if (predecessorBlock.getMarkedPosition() == predecessorBlock.getBegin()) { + partition.swapStatesAtPositions(predecessorBlock.getMarkedPosition(), predecessorPosition); + partition.swapStatesAtPositions(predecessorPosition, currentPosition + elementsToSkip + 1); + } else { + partition.swapStatesAtPositions(predecessorBlock.getMarkedPosition(), predecessorPosition); + partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.getBegin()); + partition.swapStatesAtPositions(predecessorPosition, currentPosition + elementsToSkip + 1); + } + + ++elementsToSkip; + predecessorBlock.incrementMarkedPosition(); + predecessorBlock.incrementBegin(); + } + } else { + partition.swapStates(predecessor, partition.getState(predecessorBlock.getBegin())); + predecessorBlock.incrementBegin(); + } + partition.setValue(predecessor, predecessorEntry.getValue()); + } else { + // Otherwise, we just need to update the probability for this predecessor. + partition.increaseValue(predecessor, predecessorEntry.getValue()); + } + + if (!predecessorBlock.isMarkedAsPredecessor()) { + predecessorBlocks.emplace_back(&predecessorBlock); + predecessorBlock.markAsPredecessorBlock(); + } + } + + // If we had to move some elements beyond the current element, we may have to skip them. + if (elementsToSkip > 0) { + stateIterator += elementsToSkip; + currentPosition += elementsToSkip; + } + } + + // Now we can traverse the list of states of the splitter whose predecessors we have not yet explored. + for (auto stateIterator = partition.getStates().begin() + splitter.getOriginalBegin(), stateIte = partition.getStates().begin() + splitter.getMarkedPosition(); stateIterator != stateIte; ++stateIterator) { + storm::storage::sparse::state_type currentState = *stateIterator; + + for (auto const& predecessorEntry : backwardTransitions.getRow(currentState)) { + storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn(); + Block& predecessorBlock = partition.getBlock(predecessor); + storm::storage::sparse::state_type predecessorPosition = partition.getPosition(predecessor); + + if (predecessorPosition >= predecessorBlock.getBegin()) { + partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.getBegin()); + predecessorBlock.incrementBegin(); + partition.setValue(predecessor, predecessorEntry.getValue()); + } else { + partition.increaseValue(predecessor, predecessorEntry.getValue()); + } + + if (!predecessorBlock.isMarkedAsPredecessor()) { + predecessorBlocks.emplace_back(&predecessorBlock); + predecessorBlock.markAsPredecessorBlock(); + } + } + } + + // Reset the marked position of the splitter to the begin. + splitter.setMarkedPosition(splitter.getBegin()); + + std::list blocksToSplit; + + // Now, we can iterate over the predecessor blocks and see whether we have to create a new block for + // predecessors of the splitter. + for (auto blockPtr : predecessorBlocks) { + Block& block = *blockPtr; + + block.unmarkAsPredecessorBlock(); + block.resetMarkedPosition(); + + // If we have moved the begin of the block to somewhere in the middle of the block, we need to split it. + if (block.getBegin() != block.getEnd()) { + Block& newBlock = partition.insertBlock(block); + if (!newBlock.isMarkedAsSplitter()) { + splitterQueue.push_back(&newBlock); + newBlock.markAsSplitter(); + } + + // Schedule the block of predecessors for refinement based on probabilities. + blocksToSplit.emplace_back(&newBlock); + } else { + // In this case, we can keep the block by setting its begin to the old value. + block.setBegin(block.getOriginalBegin()); + blocksToSplit.emplace_back(&block); + } + } + + // Finally, we walk through the blocks that have a transition to the splitter and split them using + // probabilistic information. + for (auto blockPtr : blocksToSplit) { + if (blockPtr->getNumberOfStates() <= 1) { + continue; + } + + refineBlockProbabilities(*blockPtr, partition, splitterQueue); + } + } + + template + template + typename DeterministicModelStrongBisimulationDecomposition::Partition DeterministicModelStrongBisimulationDecomposition::getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix const& backwardTransitions, std::string const& phiLabel, std::string const& psiLabel, bool bounded) { + std::pair statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiLabel == "true" ? storm::storage::BitVector(model.getNumberOfStates(), true) : model.getLabeledStates(phiLabel), model.getLabeledStates(psiLabel)); + Partition partition(model.getNumberOfStates(), statesWithProbability01.first, bounded ? model.getLabeledStates(psiLabel) : statesWithProbability01.second, phiLabel, psiLabel); + return partition; + } + + template + template + typename DeterministicModelStrongBisimulationDecomposition::Partition DeterministicModelStrongBisimulationDecomposition::getLabelBasedInitialPartition(ModelType const& model) { + Partition partition(model.getNumberOfStates()); + for (auto const& label : model.getStateLabeling().getAtomicPropositions()) { + if (label == "init") { + continue; + } + partition.splitLabel(model.getLabeledStates(label)); + } + return partition; + } + + template class DeterministicModelStrongBisimulationDecomposition; + template class DeterministicModelStrongBisimulationDecomposition; + } +} \ No newline at end of file diff --git a/src/storage/DeterministicModelStrongBisimulationDecomposition.h b/src/storage/DeterministicModelStrongBisimulationDecomposition.h new file mode 100644 index 000000000..e5ac547db --- /dev/null +++ b/src/storage/DeterministicModelStrongBisimulationDecomposition.h @@ -0,0 +1,452 @@ +#ifndef STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_ +#define STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_ + +#include +#include + +#include "src/storage/sparse/StateType.h" +#include "src/storage/Decomposition.h" +#include "src/models/Dtmc.h" +#include "src/models/Ctmc.h" +#include "src/storage/Distribution.h" +#include "src/utility/ConstantsComparator.h" +#include "src/utility/OsDetection.h" + +namespace storm { + namespace storage { + + /*! + * This class represents the decomposition model into its (strong) bisimulation quotient. + */ + template + class DeterministicModelStrongBisimulationDecomposition : public Decomposition { + public: + /*! + * Decomposes the given DTMC into equivalence classes under strong bisimulation. + * + * @param model The model to decompose. + * @param buildQuotient Sets whether or not the quotient model is to be built. + */ + DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc const& model, bool buildQuotient = false); + + /*! + * Decomposes the given CTMC into equivalence classes under strong bisimulation. + * + * @param model The model to decompose. + * @param buildQuotient Sets whether or not the quotient model is to be built. + */ + DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc const& model, bool buildQuotient = false); + + /*! + * Decomposes the given DTMC into equivalence classes under strong bisimulation in a way that onle safely + * preserves formulas of the form phi until psi. + * + * @param model The model to decompose. + * @param phiLabel The label that all phi states carry in the model. + * @param psiLabel The label that all psi states carry in the model. + * @param bounded If set to true, also bounded until formulas are preserved. + * @param buildQuotient Sets whether or not the quotient model is to be built. + */ + DeterministicModelStrongBisimulationDecomposition(storm::models::Dtmc const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient = false); + + /*! + * Decomposes the given CTMC into equivalence classes under strong bisimulation in a way that onle safely + * preserves formulas of the form phi until psi. + * + * @param model The model to decompose. + * @param phiLabel The label that all phi states carry in the model. + * @param psiLabel The label that all psi states carry in the model. + * @param bounded If set to true, also bounded until formulas are preserved. + * @param buildQuotient Sets whether or not the quotient model is to be built. + */ + DeterministicModelStrongBisimulationDecomposition(storm::models::Ctmc const& model, std::string const& phiLabel, std::string const& psiLabel, bool bounded, bool buildQuotient = false); + + /*! + * Retrieves the quotient of the model under the previously computed bisimulation. + * + * @return The quotient model. + */ + std::shared_ptr> getQuotient() const; + + private: + class Partition; + + class Block { + public: + typedef typename std::list::const_iterator const_iterator; + + Block(storm::storage::sparse::state_type begin, storm::storage::sparse::state_type end, Block* prev, Block* next, std::shared_ptr const& label = nullptr); + + // Prints the block. + void print(Partition const& partition) const; + + // Sets the beginning index of the block. + void setBegin(storm::storage::sparse::state_type begin); + + // Moves the beginning index of the block one step further. + void incrementBegin(); + + // Sets the end index of the block. + void setEnd(storm::storage::sparse::state_type end); + + // Moves the end index of the block one step to the front. + void decrementEnd(); + + // Returns the beginning index of the block. + storm::storage::sparse::state_type getBegin() const; + + // Returns the beginning index of the block. + storm::storage::sparse::state_type getEnd() const; + + // Retrieves the original beginning index of the block in case the begin index has been moved. + storm::storage::sparse::state_type getOriginalBegin() const; + + // Returns the iterator the block in the list of overall blocks. + const_iterator getIterator() const; + + // Returns the iterator the block in the list of overall blocks. + void setIterator(const_iterator it); + + // Returns the iterator the next block in the list of overall blocks if it exists. + const_iterator getNextIterator() const; + + // Returns the iterator the next block in the list of overall blocks if it exists. + const_iterator getPreviousIterator() const; + + // Gets the next block (if there is one). + Block& getNextBlock(); + + // Gets the next block (if there is one). + Block const& getNextBlock() const; + + // Gets a pointer to the next block (if there is one). + Block* getNextBlockPointer(); + + // Retrieves whether the block as a successor block. + bool hasNextBlock() const; + + // Gets the previous block (if there is one). + Block& getPreviousBlock(); + + // Gets a pointer to the previous block (if there is one). + Block* getPreviousBlockPointer(); + + // Gets the next block (if there is one). + Block const& getPreviousBlock() const; + + // Retrieves whether the block as a successor block. + bool hasPreviousBlock() const; + + // Checks consistency of the information in the block. + bool check() const; + + // Retrieves the number of states in this block. + std::size_t getNumberOfStates() const; + + // Checks whether the block is marked as a splitter. + bool isMarkedAsSplitter() const; + + // Marks the block as being a splitter. + void markAsSplitter(); + + // Removes the mark. + void unmarkAsSplitter(); + + // Retrieves the ID of the block. + std::size_t getId() const; + + // Retrieves the marked position in the block. + storm::storage::sparse::state_type getMarkedPosition() const; + + // Sets the marked position to the given value.. + void setMarkedPosition(storm::storage::sparse::state_type position); + + // Increases the marked position by one. + void incrementMarkedPosition(); + + // Resets the marked position to the begin of the block. + void resetMarkedPosition(); + + // Retrieves whether the block is marked as a predecessor. + bool isMarkedAsPredecessor() const; + + // Marks the block as being a predecessor block. + void markAsPredecessorBlock(); + + // Removes the marking. + void unmarkAsPredecessorBlock(); + + // Sets whether or not the block is to be interpreted as absorbing. + void setAbsorbing(bool absorbing); + + // Retrieves whether the block is to be interpreted as absorbing. + bool isAbsorbing() const; + + // Retrieves whether the block has a special label. + bool hasLabel() const; + + // Retrieves the special label of the block if it has one. + std::string const& getLabel() const; + + // Retrieves a pointer to the label of the block (which is the nullptr if there is none). + std::shared_ptr const& getLabelPtr() const; + + private: + // An iterator to itself. This is needed to conveniently insert elements in the overall list of blocks + // kept by the partition. + const_iterator selfIt; + + // Pointers to the next and previous block. + Block* next; + Block* prev; + + // The begin and end indices of the block in terms of the state vector of the partition. + storm::storage::sparse::state_type begin; + storm::storage::sparse::state_type end; + + // A field that can be used for marking the block. + bool markedAsSplitter; + + // A field that can be used for marking the block as a predecessor block. + bool markedAsPredecessorBlock; + + // A position that can be used to store a certain position within the block. + storm::storage::sparse::state_type markedPosition; + + // A flag indicating whether the block is to be interpreted as absorbing or not. + bool absorbing; + + // The ID of the block. This is only used for debugging purposes. + std::size_t id; + + // The label of the block or nullptr if it has none. + std::shared_ptr label; + + // A counter for the IDs of the blocks. + static std::size_t blockId; + }; + + class Partition { + public: + friend class Block; + + /*! + * Creates a partition with one block consisting of all the states. + * + * @param numberOfStates The number of states the partition holds. + */ + Partition(std::size_t numberOfStates); + + /*! + * Creates a partition with three blocks: one with all phi states, one with all psi states and one with + * all other states. The former two blocks are marked as being absorbing, because their outgoing + * transitions shall not be taken into account for future refinement. + * + * @param numberOfStates The number of states the partition holds. + * @param prob0States The states which have probability 0 of satisfying phi until psi. + * @param prob1States The states which have probability 1 of satisfying phi until psi. + * @param otherLabel The label that is to be attached to all other states. + * @param prob1Label The label that is to be attached to all states with probability 1. + */ + Partition(std::size_t numberOfStates, storm::storage::BitVector const& prob0States, storm::storage::BitVector const& prob1States, std::string const& otherLabel, std::string const& prob1Label); + + Partition() = default; + Partition(Partition const& other) = default; + Partition& operator=(Partition const& other) = default; +#ifndef WINDOWS + Partition(Partition&& other) = default; + Partition& operator=(Partition&& other) = default; +#endif + + /*! + * Splits all blocks of the partition such that afterwards all blocks contain only states with the label + * or no labeled state at all. + */ + void splitLabel(storm::storage::BitVector const& statesWithLabel); + + // Retrieves the size of the partition, i.e. the number of blocks. + std::size_t size() const; + + // Prints the partition to the standard output. + void print() const; + + // Splits the block at the given position and inserts a new block after the current one holding the rest + // of the states. + Block& splitBlock(Block& block, storm::storage::sparse::state_type position); + + // Inserts a block before the given block. The new block will cover all states between the beginning + // of the given block and the end of the previous block. + Block& insertBlock(Block& block); + + // Retrieves the blocks of the partition. + std::list const& getBlocks() const; + + // Retrieves the blocks of the partition. + std::list& getBlocks(); + + // Retrieves the vector of all the states. + std::vector& getStates(); + + // Checks the partition for internal consistency. + bool check() const; + + // Returns an iterator to the beginning of the states of the given block. + std::vector::iterator getBeginOfStates(Block const& block); + + // Returns an iterator to the beginning of the states of the given block. + std::vector::iterator getEndOfStates(Block const& block); + + // Returns an iterator to the beginning of the states of the given block. + std::vector::const_iterator getBeginOfStates(Block const& block) const; + + // Returns an iterator to the beginning of the states of the given block. + std::vector::const_iterator getEndOfStates(Block const& block) const; + + // Returns an iterator to the beginning of the states of the given block. + typename std::vector::iterator getBeginOfValues(Block const& block); + + // Returns an iterator to the beginning of the states of the given block. + typename std::vector::iterator getEndOfValues(Block const& block); + + // Swaps the positions of the two given states. + void swapStates(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2); + + // Swaps the positions of the two states given by their positions. + void swapStatesAtPositions(storm::storage::sparse::state_type position1, storm::storage::sparse::state_type position2); + + // Retrieves the block of the given state. + Block& getBlock(storm::storage::sparse::state_type state); + + // Retrieves the block of the given state. + Block const& getBlock(storm::storage::sparse::state_type state) const; + + // Retrieves the position of the given state. + storm::storage::sparse::state_type const& getPosition(storm::storage::sparse::state_type state) const; + + // Retrieves the position of the given state. + void setPosition(storm::storage::sparse::state_type state, storm::storage::sparse::state_type position); + + // Sets the position of the state to the given position. + storm::storage::sparse::state_type const& getState(storm::storage::sparse::state_type position) const; + + // Retrieves the value for the given state. + ValueType const& getValue(storm::storage::sparse::state_type state) const; + + // Retrieves the value at the given position. + ValueType const& getValueAtPosition(storm::storage::sparse::state_type position) const; + + // Sets the given value for the given state. + void setValue(storm::storage::sparse::state_type state, ValueType value); + + // Retrieves the vector with the probabilities going into the current splitter. + std::vector& getValues(); + + // Increases the value for the given state by the specified amount. + void increaseValue(storm::storage::sparse::state_type state, ValueType value); + + // Updates the block mapping for the given range of states to the specified block. + void updateBlockMapping(Block& block, std::vector::iterator first, std::vector::iterator end); + + // Retrieves the first block of the partition. + Block& getFirstBlock(); + private: + // The list of blocks in the partition. + std::list blocks; + + // A mapping of states to their blocks. + std::vector stateToBlockMapping; + + // A vector containing all the states. It is ordered in a special way such that the blocks only need to + // define their start/end indices. + std::vector states; + + // This vector keeps track of the position of each state in the state vector. + std::vector positions; + + // This vector stores the probabilities of going to the current splitter. + std::vector values; + }; + + /*! + * Performs the partition refinement on the model and thereby computes the equivalence classes under strong + * bisimulation equivalence. If required, the quotient model is built and may be retrieved using + * getQuotient(). + * + * @param model The model on whose state space to compute the coarses strong bisimulation relation. + * @param backwardTransitions The backward transitions of the model. + * @param The initial partition. + * @param buildQuotient If set, the quotient model is built and may be retrieved using the getQuotient() + * method. + */ + template + void partitionRefinement(ModelType const& model, storm::storage::SparseMatrix const& backwardTransitions, Partition& partition, bool buildQuotient); + + /*! + * Refines the partition based on the provided splitter. After calling this method all blocks are stable + * with respect to the splitter. + * + * @param backwardTransitions A matrix that can be used to retrieve the predecessors (and their + * probabilities). + * @param splitter The splitter to use. + * @param partition The partition to split. + * @param splitterQueue A queue into which all blocks that were split are inserted so they can be treated + * as splitters in the future. + */ + void refinePartition(storm::storage::SparseMatrix const& backwardTransitions, Block& splitter, Partition& partition, std::deque& splitterQueue); + + /*! + * Refines the block based on their probability values (leading into the splitter). + * + * @param block The block to refine. + * @param partition The partition that contains the block. + * @param splitterQueue A queue into which all blocks that were split are inserted so they can be treated + * as splitters in the future. + */ + void refineBlockProbabilities(Block& block, Partition& partition, std::deque& splitterQueue); + + /*! + * Builds the quotient model based on the previously computed equivalence classes (stored in the blocks + * of the decomposition. + * + * @param model The model whose state space was used for computing the equivalence classes. This is used for + * determining the transitions of each equivalence class. + * @param partition The previously computed partition. This is used for quickly retrieving the block of a + * state. + */ + template + void buildQuotient(ModelType const& model, Partition const& partition); + + /*! + * Creates the measure-driven initial partition for reaching psi states from phi states. + * + * @param model The model whose state space is partitioned based on reachability of psi states from phi + * states. + * @param backwardTransitions The backward transitions of the model. + * @param phiLabel The label that all phi states carry in the model. + * @param psiLabel The label that all psi states carry in the model. + * @param bounded If set to true, the initial partition will be chosen in such a way that preserves bounded + * reachability queries. + * @return The resulting partition. + */ + template + Partition getMeasureDrivenInitialPartition(ModelType const& model, storm::storage::SparseMatrix const& backwardTransitions, std::string const& phiLabel, std::string const& psiLabel, bool bounded = false); + + /*! + * Creates the initial partition based on all the labels in the given model. + * + * @param model The model whose state space is partitioned based on its labels. + * @return The resulting partition. + */ + template + Partition getLabelBasedInitialPartition(ModelType const& model); + + // If required, a quotient model is built and stored in this member. + std::shared_ptr> quotient; + + // A comparator that is used for determining whether two probabilities are considered to be equal. + storm::utility::ConstantsComparator comparator; + }; + } +} + +#endif /* STORM_STORAGE_DETERMINISTICMODELSTRONGBISIMULATIONDECOMPOSITION_H_ */ \ No newline at end of file diff --git a/src/storage/Distribution.cpp b/src/storage/Distribution.cpp new file mode 100644 index 000000000..50e93f567 --- /dev/null +++ b/src/storage/Distribution.cpp @@ -0,0 +1,100 @@ +#include "src/storage/Distribution.h" + +#include +#include + +#include "src/settings/SettingsManager.h" + +namespace storm { + namespace storage { + + template + Distribution::Distribution() : hash(0) { + // Intentionally left empty. + } + + template + bool Distribution::operator==(Distribution const& other) const { + // We need to check equality by ourselves, because we need to account for epsilon differences. + if (this->distribution.size() != other.distribution.size() || this->getHash() != other.getHash()) { + return false; + } + + auto first1 = this->distribution.begin(); + auto last1 = this->distribution.end(); + auto first2 = other.distribution.begin(); + auto last2 = other.distribution.end(); + + for (; first1 != last1; ++first1, ++first2) { + if (first1->first != first2->first) { + return false; + } + if (std::abs(first1->second - first2->second) > 1e-6) { + return false; + } + } + + return true; + } + + template + void Distribution::addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability) { + if (this->distribution.find(state) == this->distribution.end()) { + this->hash += static_cast(state); + } + this->distribution[state] += probability; + } + + template + typename Distribution::iterator Distribution::begin() { + return this->distribution.begin(); + } + + template + typename Distribution::const_iterator Distribution::begin() const { + return this->distribution.begin(); + } + + template + typename Distribution::iterator Distribution::end() { + return this->distribution.end(); + } + + template + typename Distribution::const_iterator Distribution::end() const { + return this->distribution.end(); + } + + template + void Distribution::scale(storm::storage::sparse::state_type const& state) { + auto probabilityIterator = this->distribution.find(state); + if (probabilityIterator != this->distribution.end()) { + ValueType scaleValue = 1 / probabilityIterator->second; + this->distribution.erase(probabilityIterator); + + for (auto& entry : this->distribution) { + entry.second *= scaleValue; + } + } + } + + template + std::size_t Distribution::getHash() const { + return this->hash ^ (this->distribution.size() << 8); + } + + template + std::ostream& operator<<(std::ostream& out, Distribution const& distribution) { + out << "{"; + for (auto const& entry : distribution) { + out << "[" << entry.second << ": " << entry.first << "], "; + } + out << "}"; + + return out; + } + + template class Distribution; + template std::ostream& operator<<(std::ostream& out, Distribution const& distribution); + } +} diff --git a/src/storage/Distribution.h b/src/storage/Distribution.h new file mode 100644 index 000000000..010cf43bc --- /dev/null +++ b/src/storage/Distribution.h @@ -0,0 +1,110 @@ +#ifndef STORM_STORAGE_DISTRIBUTION_H_ +#define STORM_STORAGE_DISTRIBUTION_H_ + +#include +#include +#include + +#include "src/storage/sparse/StateType.h" + +namespace storm { + namespace storage { + + template + class Distribution { + public: + typedef boost::container::flat_map container_type; + typedef typename container_type::iterator iterator; + typedef typename container_type::const_iterator const_iterator; + + /*! + * Creates an empty distribution. + */ + Distribution(); + + /*! + * Checks whether the two distributions specify the same probabilities to go to the same states. + * + * @param other The distribution with which the current distribution is to be compared. + * @return True iff the two distributions are equal. + */ + bool operator==(Distribution const& other) const; + + /*! + * Assigns the given state the given probability under this distribution. + * + * @param state The state to which to assign the probability. + * @param probability The probability to assign. + */ + void addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability); + + /*! + * Retrieves an iterator to the elements in this distribution. + * + * @return The iterator to the elements in this distribution. + */ + iterator begin(); + + /*! + * Retrieves an iterator to the elements in this distribution. + * + * @return The iterator to the elements in this distribution. + */ + const_iterator begin() const; + + /*! + * Retrieves an iterator past the elements in this distribution. + * + * @return The iterator past the elements in this distribution. + */ + iterator end(); + + /*! + * Retrieves an iterator past the elements in this distribution. + * + * @return The iterator past the elements in this distribution. + */ + const_iterator end() const; + + /*! + * Scales the distribution by multiplying all the probabilities with 1/p where p is the probability of moving + * to the given state and sets the probability of moving to the given state to zero. If the probability is + * already zero, this operation has no effect. + * + * @param state The state whose associated probability is used to scale the distribution. + */ + void scale(storm::storage::sparse::state_type const& state); + + /*! + * Retrieves the hash value of the distribution. + * + * @return The hash value of the distribution. + */ + std::size_t getHash() const; + + private: + // A list of states and the probabilities that are assigned to them. + container_type distribution; + + // A hash value that is maintained to allow for quicker equality comparison between distribution.s + std::size_t hash; + }; + + template + std::ostream& operator<<(std::ostream& out, Distribution const& distribution); + } +} + +namespace std { + + template + struct hash> { + std::size_t operator()(storm::storage::Distribution const& distribution) const { + return (distribution.getHash()); + } + }; + +} + + +#endif /* STORM_STORAGE_DISTRIBUTION_H_ */ \ No newline at end of file diff --git a/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp b/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp new file mode 100644 index 000000000..bf41414de --- /dev/null +++ b/src/storage/NaiveDeterministicModelBisimulationDecomposition.cpp @@ -0,0 +1,262 @@ +#include "src/storage/NaiveDeterministicModelBisimulationDecomposition.h" + +#include +#include + +namespace storm { + namespace storage { + + template + NaiveDeterministicModelBisimulationDecomposition::NaiveDeterministicModelBisimulationDecomposition(storm::models::Dtmc const& dtmc, bool weak) { + computeBisimulationEquivalenceClasses(dtmc, weak); + } + + template + void NaiveDeterministicModelBisimulationDecomposition::computeBisimulationEquivalenceClasses(storm::models::Dtmc const& dtmc, bool weak) { + std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now(); + // We start by computing the initial partition. In particular, we also keep a mapping of states to their blocks. + std::vector stateToBlockMapping(dtmc.getNumberOfStates()); + storm::storage::BitVector labeledStates = dtmc.getLabeledStates("observe0Greater1"); + this->blocks.emplace_back(labeledStates.begin(), labeledStates.end()); + std::for_each(labeledStates.begin(), labeledStates.end(), [&] (storm::storage::sparse::state_type const& state) { stateToBlockMapping[state] = 0; } ); + labeledStates.complement(); + this->blocks.emplace_back(labeledStates.begin(), labeledStates.end()); + std::for_each(labeledStates.begin(), labeledStates.end(), [&] (storm::storage::sparse::state_type const& state) { stateToBlockMapping[state] = 1; } ); + + // Check whether any of the blocks is empty and remove it. + auto eraseIterator = std::remove_if(this->blocks.begin(), this->blocks.end(), [] (typename NaiveDeterministicModelBisimulationDecomposition::block_type const& a) { return a.empty(); }); + this->blocks.erase(eraseIterator, this->blocks.end()); + + // Create empty distributions for the two initial blocks. + std::vector> distributions(2); + + // Retrieve the backward transitions to allow for better checking of states that need to be re-examined. + storm::storage::SparseMatrix const& backwardTransitions = dtmc.getBackwardTransitions(); + + // Initially, both blocks are potential splitters. A splitter is marked as a pair in which the first entry + // is the ID of the parent block of the splitter and the second entry is the block ID of the splitter itself. + std::deque refinementQueue; + storm::storage::BitVector blocksInRefinementQueue(this->size()); + for (uint_fast64_t index = 0; index < this->blocks.size(); ++index) { + refinementQueue.push_back(index); + } + + // As long as there are blocks to refine, well, refine them. + uint_fast64_t iteration = 0; + while (!refinementQueue.empty()) { + ++iteration; + + // Optionally sort the queue to potentially speed up the convergence. + // std::sort(refinementQueue.begin(), refinementQueue.end(), [=] (std::size_t const& a, std::size_t const& b) { return this->blocks[a].size() > this->blocks[b].size(); }); + + std::size_t currentBlock = refinementQueue.front(); + refinementQueue.pop_front(); + blocksInRefinementQueue.set(currentBlock, false); + + splitBlock(dtmc, backwardTransitions, currentBlock, stateToBlockMapping, distributions, blocksInRefinementQueue, refinementQueue, weak); + } + + std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart; + + std::cout << "Bisimulation quotient has " << this->blocks.size() << " blocks and took " << iteration << " iterations and " << std::chrono::duration_cast(totalTime).count() << "ms." << std::endl; + } + + template + std::size_t NaiveDeterministicModelBisimulationDecomposition::splitPredecessorsGraphBased(storm::models::Dtmc const& dtmc, storm::storage::SparseMatrix const& backwardTransitions, std::size_t const& blockId, std::vector& stateToBlockMapping, std::vector>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque& graphRefinementQueue, storm::storage::BitVector& splitBlocks) { +// std::cout << "entering " << std::endl; + std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now(); + +// std::cout << "refining predecessors of " << blockId << std::endl; + + // This method tries to split the blocks of predecessors of the provided block by checking whether there is + // a transition into the current block or not. + std::unordered_map::block_type> predecessorBlockToNewBlock; + + // Now for each predecessor block which state could actually reach the current block. + std::chrono::high_resolution_clock::time_point predIterStart = std::chrono::high_resolution_clock::now(); + int predCounter = 0; + storm::storage::BitVector preds(dtmc.getNumberOfStates()); + for (auto const& state : this->blocks[blockId]) { + for (auto const& predecessorEntry : backwardTransitions.getRow(state)) { + ++predCounter; + storm::storage::sparse::state_type predecessor = predecessorEntry.getColumn(); + if (!preds.get(predecessor)) { +// std::cout << "having pred " << predecessor << " with block " << stateToBlockMapping[predecessor] << std::endl; + predecessorBlockToNewBlock[stateToBlockMapping[predecessor]].insert(predecessor); + preds.set(predecessor); + } + } + } + std::chrono::high_resolution_clock::duration predIterTime = std::chrono::high_resolution_clock::now() - predIterStart; + std::cout << "took " << std::chrono::duration_cast(predIterTime).count() << "ms. for " << predCounter << " predecessors" << std::endl; + + // Now, we can check for each predecessor block whether it needs to be split. + for (auto const& blockNewBlockPair : predecessorBlockToNewBlock) { + if (this->blocks[blockNewBlockPair.first].size() > blockNewBlockPair.second.size()) { + std::cout << "splitting!" << std::endl; +// std::cout << "found that not all " << this->blocks[blockNewBlockPair.first].size() << " states of block " << blockNewBlockPair.first << " have a successor in " << blockId << " but only " << blockNewBlockPair.second.size() << std::endl; +// std::cout << "original: " << this->blocks[blockNewBlockPair.first] << std::endl; +// std::cout << "states with pred: " << blockNewBlockPair.second << std::endl; + // Now update the block mapping for the smaller of the two blocks. + typename NaiveDeterministicModelBisimulationDecomposition::block_type smallerBlock; + typename NaiveDeterministicModelBisimulationDecomposition::block_type largerBlock; + if (blockNewBlockPair.second.size() < this->blocks[blockNewBlockPair.first].size()/2) { + smallerBlock = std::move(blockNewBlockPair.second); + std::set_difference(this->blocks[blockNewBlockPair.first].begin(), this->blocks[blockNewBlockPair.first].end(), smallerBlock.begin(), smallerBlock.end(), std::inserter(largerBlock, largerBlock.begin())); + } else { + largerBlock = std::move(blockNewBlockPair.second); + std::set_difference(this->blocks[blockNewBlockPair.first].begin(), this->blocks[blockNewBlockPair.first].end(), largerBlock.begin(), largerBlock.end(), std::inserter(smallerBlock, smallerBlock.begin())); + } + +// std::cout << "created a smaller block with " << smallerBlock.size() << " states and a larger one with " << largerBlock.size() << "states" << std::endl; +// std::cout << "smaller: " << smallerBlock << std::endl; +// std::cout << "larger: " << largerBlock << std::endl; + + this->blocks.emplace_back(std::move(smallerBlock)); + this->blocks[blockNewBlockPair.first] = std::move(largerBlock); + + // Update the block mapping of all moved states. + std::size_t newBlockId = this->blocks.size() - 1; + for (auto const& state : this->blocks.back()) { + stateToBlockMapping[state] = newBlockId; +// std::cout << "updating " << state << " to block " << newBlockId << std::endl; + } + + blocksInRefinementQueue.resize(this->blocks.size()); + + // Add both parts of the old block to the queue. + if (!blocksInRefinementQueue.get(blockNewBlockPair.first)) { + graphRefinementQueue.push_back(blockNewBlockPair.first); + blocksInRefinementQueue.set(blockNewBlockPair.first, true); +// std::cout << "adding " << blockNewBlockPair.first << " to refine further using graph-based analysis " << std::endl; + } + + graphRefinementQueue.push_back(this->blocks.size() - 1); + blocksInRefinementQueue.set(this->blocks.size() - 1); +// std::cout << "adding " << this->blocks.size() - 1 << " to refine further using graph-based analysis " << std::endl; + } + } + std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart; + std::cout << "refinement of predecessors of block " << blockId << " took " << std::chrono::duration_cast(totalTime).count() << "ms." << std::endl; + +// std::cout << "done. "<< std::endl; + + // FIXME + return 0; + } + + template + std::size_t NaiveDeterministicModelBisimulationDecomposition::splitBlock(storm::models::Dtmc const& dtmc, storm::storage::SparseMatrix const& backwardTransitions, std::size_t const& blockId, std::vector& stateToBlockMapping, std::vector>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque& refinementQueue, bool weakBisimulation) { + std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now(); + std::unordered_map, typename NaiveDeterministicModelBisimulationDecomposition::block_type> distributionToNewBlocks; + + // Traverse all states of the block and check whether they have different distributions. + std::chrono::high_resolution_clock::time_point gatherStart = std::chrono::high_resolution_clock::now(); + for (auto const& state : this->blocks[blockId]) { + // Now construct the distribution of this state wrt. to the current partitioning. + storm::storage::Distribution distribution; + for (auto const& successorEntry : dtmc.getTransitionMatrix().getRow(state)) { + distribution.addProbability(static_cast(stateToBlockMapping[successorEntry.getColumn()]), successorEntry.getValue()); + } + + // If we are requested to compute a weak bisimulation, we need to scale the distribution with the + // self-loop probability. + if (weakBisimulation) { + distribution.scale(blockId); + } + + // If the distribution already exists, we simply add the state. Otherwise, we open a new block. + auto distributionIterator = distributionToNewBlocks.find(distribution); + if (distributionIterator != distributionToNewBlocks.end()) { + distributionIterator->second.insert(state); + } else { + distributionToNewBlocks[distribution].insert(state); + } + } + + std::chrono::high_resolution_clock::duration gatherTime = std::chrono::high_resolution_clock::now() - gatherStart; + std::cout << "time to iterate over all states was " << std::chrono::duration_cast(gatherTime).count() << "ms." << std::endl; + + // Now we are ready to split the block. + if (distributionToNewBlocks.size() == 1) { + // If there is just one behavior, we just set the distribution as the new one for this block. + // distributions[blockId] = std::move(distributionToNewBlocks.begin()->first); + } else { + // In this case, we need to split the block. + typename NaiveDeterministicModelBisimulationDecomposition::block_type tmpBlock; + + auto distributionIterator = distributionToNewBlocks.begin(); + STORM_LOG_ASSERT(distributionIterator != distributionToNewBlocks.end(), "One block has no distribution..."); + + // distributions[blockId] = std::move(distributionIterator->first); + tmpBlock = std::move(distributionIterator->second); + std::swap(this->blocks[blockId], tmpBlock); + ++distributionIterator; + + // Remember the number of blocks prior to splitting for later use. + std::size_t beforeNumberOfBlocks = this->blocks.size(); + + for (; distributionIterator != distributionToNewBlocks.end(); ++distributionIterator) { + // In this case, we need to move the newly created block to the end of the list of actual blocks. + this->blocks.emplace_back(std::move(distributionIterator->second)); + distributions.emplace_back(std::move(distributionIterator->first)); + + // Update the mapping of states to their blocks. + std::size_t newBlockId = this->blocks.size() - 1; + for (auto const& state : this->blocks.back()) { + stateToBlockMapping[state] = newBlockId; + } + } + + // Now that we have split the block, we try to trigger a chain of graph-based splittings. That is, we + // try to split the predecessors of the current block by checking whether they have some transition + // to one given sub-block of the current block. + std::deque localRefinementQueue; + storm::storage::BitVector blocksInLocalRefinementQueue(this->size()); + localRefinementQueue.push_back(blockId); +// std::cout << "adding " << blockId << " to local ref queue " << std::endl; + blocksInLocalRefinementQueue.set(blockId); + for (std::size_t i = beforeNumberOfBlocks; i < this->blocks.size(); ++i) { + localRefinementQueue.push_back(i); + blocksInLocalRefinementQueue.set(i); +// std::cout << "adding " << i << " to local refinement queue " << std::endl; + } + + // We need to keep track of which blocks were split so we can later add all their predecessors as + // candidates for furter splitting. + storm::storage::BitVector splitBlocks(this->size()); + +// while (!localRefinementQueue.empty()) { +// std::size_t currentBlock = localRefinementQueue.front(); +// localRefinementQueue.pop_front(); +// blocksInLocalRefinementQueue.set(currentBlock, false); +// +// splitPredecessorsGraphBased(dtmc, backwardTransitions, currentBlock, stateToBlockMapping, distributions, blocksInLocalRefinementQueue, localRefinementQueue, splitBlocks); +// } + +// std::cout << "split blocks: " << std::endl; +// std::cout << splitBlocks << std::endl; + + // Since we created some new blocks, we need to extend the bit vector storing the blocks in the + // refinement queue. + blocksInRefinementQueue.resize(this->blocks.size()); + + // Insert blocks that possibly need a refinement into the queue. + for (auto const& state : tmpBlock) { + for (auto const& predecessor : backwardTransitions.getRow(state)) { + if (!blocksInRefinementQueue.get(stateToBlockMapping[predecessor.getColumn()])) { + blocksInRefinementQueue.set(stateToBlockMapping[predecessor.getColumn()]); + refinementQueue.push_back(stateToBlockMapping[predecessor.getColumn()]); + } + } + } + } + + std::chrono::high_resolution_clock::duration totalTime = std::chrono::high_resolution_clock::now() - totalStart; + std::cout << "refinement of block " << blockId << " took " << std::chrono::duration_cast(totalTime).count() << "ms." << std::endl; + return distributionToNewBlocks.size(); + } + + template class NaiveDeterministicModelBisimulationDecomposition; + } +} \ No newline at end of file diff --git a/src/storage/NaiveDeterministicModelBisimulationDecomposition.h b/src/storage/NaiveDeterministicModelBisimulationDecomposition.h new file mode 100644 index 000000000..300922cf5 --- /dev/null +++ b/src/storage/NaiveDeterministicModelBisimulationDecomposition.h @@ -0,0 +1,35 @@ +#ifndef STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ +#define STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ + +#include +#include + +#include "src/storage/Decomposition.h" +#include "src/models/Dtmc.h" +#include "src/storage/Distribution.h" + +namespace storm { + namespace storage { + + /*! + * This class represents the decomposition model into its bisimulation quotient. + */ + template + class NaiveDeterministicModelBisimulationDecomposition : public Decomposition { + public: + NaiveDeterministicModelBisimulationDecomposition() = default; + + /*! + * Decomposes the given DTMC into equivalence classes under weak or strong bisimulation. + */ + NaiveDeterministicModelBisimulationDecomposition(storm::models::Dtmc const& model, bool weak = false); + + private: + void computeBisimulationEquivalenceClasses(storm::models::Dtmc const& model, bool weak); + std::size_t splitPredecessorsGraphBased(storm::models::Dtmc const& dtmc, storm::storage::SparseMatrix const& backwardTransitions, std::size_t const& block, std::vector& stateToBlockMapping, std::vector>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque& refinementQueue, storm::storage::BitVector& splitBlocks); + std::size_t splitBlock(storm::models::Dtmc const& dtmc, storm::storage::SparseMatrix const& backwardTransitions, std::size_t const& block, std::vector& stateToBlockMapping, std::vector>& distributions, storm::storage::BitVector& blocksInRefinementQueue, std::deque& refinementQueue, bool weakBisimulation); + }; + } +} + +#endif /* STORM_STORAGE_NAIVEDETERMINISTICMODELBISIMULATIONDECOMPOSITION_H_ */ \ No newline at end of file diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 6d9531428..0b84e9936 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -239,7 +239,7 @@ namespace storm { template SparseMatrix::SparseMatrix(index_type columnCount, std::vector const& rowIndications, std::vector> const& columnsAndValues, std::vector const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) { for (auto const& element : *this) { - if (element.getValue() != storm::utility::constantZero()) { + if (!storm::utility::isZero(element.getValue())) { ++this->nonzeroEntryCount; } } @@ -248,7 +248,7 @@ namespace storm { template SparseMatrix::SparseMatrix(index_type columnCount, std::vector&& rowIndications, std::vector>&& columnsAndValues, std::vector&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) { for (auto const& element : *this) { - if (element.getValue() != storm::utility::constantZero()) { + if (!storm::utility::isZero(element.getValue())) { ++this->nonzeroEntryCount; } } @@ -604,7 +604,7 @@ namespace storm { // First, we need to count how many entries each column has. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::constantZero()) { + if (!storm::utility::isZero(transition.getValue())) { ++rowIndications[transition.getColumn() + 1]; } } @@ -623,7 +623,7 @@ namespace storm { // Now we are ready to actually fill in the values of the transposed matrix. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::constantZero()) { + if (!storm::utility::isZero(transition.getValue())) { columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); nextIndices[transition.getColumn()]++; } diff --git a/src/storage/StateBlock.cpp b/src/storage/StateBlock.cpp index fb4f21133..bc270872d 100644 --- a/src/storage/StateBlock.cpp +++ b/src/storage/StateBlock.cpp @@ -30,6 +30,10 @@ namespace storm { states.insert(state); } + StateBlock::iterator StateBlock::insert(container_type::const_iterator iterator, value_type const& state) { + return states.insert(iterator, state); + } + void StateBlock::erase(value_type const& state) { states.erase(state); } @@ -55,5 +59,6 @@ namespace storm { out << block.getStates(); return out; } + } } \ No newline at end of file diff --git a/src/storage/StateBlock.h b/src/storage/StateBlock.h index dcd516277..e9a8d5ebe 100644 --- a/src/storage/StateBlock.h +++ b/src/storage/StateBlock.h @@ -2,6 +2,7 @@ #define STORM_STORAGE_BLOCK_H_ #include +#include #include "src/utility/OsDetection.h" #include "src/storage/sparse/StateType.h" @@ -35,10 +36,15 @@ namespace storm { * * @param first The first element of the range to insert. * @param last The last element of the range (that is itself not inserted). + * @param sortedAndUnique If set to true, the input range is assumed to be sorted and duplicate-free. */ template - StateBlock(InputIterator first, InputIterator last) : states(first, last) { - // Intentionally left empty. + StateBlock(InputIterator first, InputIterator last, bool sortedAndUnique = false) { + if (sortedAndUnique) { + this->states = container_type(boost::container::ordered_unique_range_t(), first, last); + } else { + this->states = container_type(first, last); + } } /*! @@ -101,6 +107,13 @@ namespace storm { * @param state The state to add to this SCC. */ void insert(value_type const& state); + + /*! + * Inserts the given element into this SCC. + * + * @param state The state to add to this SCC. + */ + iterator insert(container_type::const_iterator iterator, value_type const& state); /*! * Removes the given element from this SCC. diff --git a/src/storage/StronglyConnectedComponentDecomposition.cpp b/src/storage/StronglyConnectedComponentDecomposition.cpp index 7c6988303..63b6d1de5 100644 --- a/src/storage/StronglyConnectedComponentDecomposition.cpp +++ b/src/storage/StronglyConnectedComponentDecomposition.cpp @@ -171,7 +171,7 @@ namespace storm { p.push_back(currentState); for (auto const& successor : transitionMatrix.getRowGroup(currentState)) { - if (subsystem.get(successor.getColumn()) && successor.getValue() != storm::utility::constantZero()) { + if (subsystem.get(successor.getColumn()) && !storm::utility::isZero(successor.getValue())) { if (currentState == successor.getColumn()) { statesWithSelfLoop.set(currentState); } diff --git a/src/storage/expressions/ExpressionEvaluation.h b/src/storage/expressions/ExpressionEvaluation.h index 7b35109f7..89a32b038 100644 --- a/src/storage/expressions/ExpressionEvaluation.h +++ b/src/storage/expressions/ExpressionEvaluation.h @@ -1,4 +1,4 @@ -/** +/** * @file: ExpressionEvaluation.h * @author: Sebastian Junges * @@ -20,56 +20,65 @@ #include "src/storage/parameters.h" namespace storm { -namespace expressions { - - template - struct StateType - { - typedef int type; - }; - + namespace expressions { + + template + struct StateType + { + typedef int type; + }; + #ifdef PARAMETRIC_SYSTEMS - template<> - struct StateType - { - typedef std::map type; - }; - - template<> - struct StateType - { - typedef std::map type; - }; + template<> + struct StateType + { + typedef std::map type; + }; + + template<> + struct StateType + { + typedef std::map type; + }; #endif - template - class ExpressionEvaluationVisitor : public ExpressionVisitor - { + template + class ExpressionEvaluationVisitor : public ExpressionVisitor + { public: ExpressionEvaluationVisitor(S* sharedState) - : mSharedState(sharedState) + : mSharedState(sharedState), cache(new carl::Cache>()) { - + } + + ExpressionEvaluationVisitor(S* sharedState, std::shared_ptr>> cache) + : mSharedState(sharedState), cache(cache) + { + + } + + + virtual ~ExpressionEvaluationVisitor() { + + } - virtual ~ExpressionEvaluationVisitor() {} - - virtual void visit(IfThenElseExpression const* expression) + virtual void visit(IfThenElseExpression const* expression) { std::cout << "ite" << std::endl; } - - virtual void visit(BinaryBooleanFunctionExpression const* expression) + + virtual void visit(BinaryBooleanFunctionExpression const* expression) { std::cout << "bbf" << std::endl; } - virtual void visit(BinaryNumericalFunctionExpression const* expression) + virtual void visit(BinaryNumericalFunctionExpression const* expression) { - ExpressionEvaluationVisitor* visitor = new ExpressionEvaluationVisitor(mSharedState); + ExpressionEvaluationVisitor* visitor = new ExpressionEvaluationVisitor(mSharedState, this->cache); expression->getFirstOperand()->accept(visitor); mValue = visitor->value(); - expression->getSecondOperand()->accept(visitor); + expression->getSecondOperand()->accept(visitor); switch(expression->getOperatorType()) { case BinaryNumericalFunctionExpression::OperatorType::Plus: @@ -91,102 +100,110 @@ namespace expressions { delete visitor; } - virtual void visit(BinaryRelationExpression const* expression) + virtual void visit(BinaryRelationExpression const* expression) { std::cout << "br" << std::endl; } - virtual void visit(VariableExpression const* expression) + virtual void visit(VariableExpression const* expression) { - std::string const& varName= expression->getVariableName(); + std::string const& varName= expression->getVariableName(); auto it = mSharedState->find(varName); if(it != mSharedState->end()) { - mValue = T(it->second); + mValue = convertVariableToPolynomial(it->second); +// mValue = T(typename T::PolyType(typename T::PolyType::PolyType(it->second), cache)); } else { carl::Variable nVar = carl::VariablePool::getInstance().getFreshVariable(varName); mSharedState->emplace(varName,nVar); - mValue = T(nVar); + mValue = convertVariableToPolynomial(nVar); } } - virtual void visit(UnaryBooleanFunctionExpression const* expression) + virtual void visit(UnaryBooleanFunctionExpression const* expression) { std::cout << "ubf" << std::endl; } - virtual void visit(UnaryNumericalFunctionExpression const* expression) + virtual void visit(UnaryNumericalFunctionExpression const* expression) { std::cout << "unf" << std::endl; } - virtual void visit(BooleanLiteralExpression const* expression) + virtual void visit(BooleanLiteralExpression const* expression) { std::cout << "bl" << std::endl; } - virtual void visit(IntegerLiteralExpression const* expression) + virtual void visit(IntegerLiteralExpression const* expression) { - // mValue = T(typename T::CoeffType(std::to_string(expression->getValue()), 10)); - mValue = T(expression->getValue()); + mValue = T(typename T::PolyType(typename T::CoeffType(expression->getValue()))); } - virtual void visit(DoubleLiteralExpression const* expression) + virtual void visit(DoubleLiteralExpression const* expression) { std::stringstream str; str << std::fixed << std::setprecision( 3 ) << expression->getValue(); - mValue = T(carl::rationalize(str.str())); - + mValue = T(carl::rationalize(str.str())); } - - const T& value() const - { - return mValue; - } - + + template> = carl::dummy> + T convertVariableToPolynomial(carl::Variable const& nVar) { + return T(typename T::PolyType(typename T::PolyType::PolyType(nVar), cache)); + } + + template> = carl::dummy> + T convertVariableToPolynomial(carl::Variable const& nVar) { + return T(nVar); + } + + const T& value() const + { + return mValue; + } + private: - S* mSharedState; - T mValue; - }; - - template - class ExpressionEvaluation - { + S* mSharedState; + T mValue; + std::shared_ptr>> cache; + }; + + template + class ExpressionEvaluation + { public: - ExpressionEvaluation() : mState() - { - - } - - - T evaluate(Expression const& expr, storm::expressions::SimpleValuation const* val) - { - ExpressionEvaluationVisitor::type>* visitor = new ExpressionEvaluationVisitor::type>(&mState); - Expression expressionToTranslate = expr.substitute(*val); - expressionToTranslate.getBaseExpression().accept(visitor); - T result = visitor->value(); -// result.simplify(); - delete visitor; - return result; - } - + ExpressionEvaluation() : mState(), cache(new carl::Cache>(100000)) + { + // Intentionally left empty. + } + + + T evaluate(Expression const& expr, storm::expressions::SimpleValuation const* val) + { + ExpressionEvaluationVisitor::type>* visitor = new ExpressionEvaluationVisitor::type>(&mState, cache); + Expression expressionToTranslate = expr.substitute(*val); + expressionToTranslate.getBaseExpression().accept(visitor); + T result = visitor->value(); + // result.simplify(); + delete visitor; + return result; + } + protected: typename StateType::type mState; - }; - - /** - * For doubles, we keep using the getValueAs from the expressions, as this should be more efficient. - */ - template<> - class ExpressionEvaluation - { + std::shared_ptr>> cache; + }; + + /** + * For doubles, we keep using the getValueAs from the expressions, as this should be more efficient. + */ + template<> + class ExpressionEvaluation + { public: - double evaluate(Expression const& expr, storm::expressions::SimpleValuation const* val) const - { - return expr.evaluateAsDouble(val); - } - }; - - - - -} + double evaluate(Expression const& expr, storm::expressions::SimpleValuation const* val) const + { + return expr.evaluateAsDouble(val); + } + }; + + } } #endif \ No newline at end of file diff --git a/src/storage/parameters.h b/src/storage/parameters.h index 923ef5fc9..5313a0e23 100644 --- a/src/storage/parameters.h +++ b/src/storage/parameters.h @@ -10,9 +10,12 @@ namespace storm { - typedef carl::MultivariatePolynomial Polynomial; +// typedef carl::MultivariatePolynomial Polynomial; typedef carl::Variable Variable; - + typedef carl::MultivariatePolynomial RawPolynomial; +// typedef carl::FactorizedPolynomial Polynomial; + + typedef RawPolynomial Polynomial; typedef carl::CompareRelation CompareRelation; typedef carl::RationalFunction RationalFunction; diff --git a/src/stormParametric.cpp b/src/stormParametric.cpp index f5f0271a0..109cebb8a 100644 --- a/src/stormParametric.cpp +++ b/src/stormParametric.cpp @@ -1,6 +1,7 @@ // Include generated headers. #include "storm-config.h" -#include "storm-version.h" + +#include // Include other headers. #include "src/exceptions/BaseException.h" @@ -9,14 +10,124 @@ #include "src/utility/export.h" #include "src/modelchecker/reachability/CollectConstraints.h" -#include "src/modelchecker/reachability/DirectEncoding.h" +//#include "src/modelchecker/reachability/DirectEncoding.h" +#include "src/storage/DeterministicModelStrongBisimulationDecomposition.h" #include "src/modelchecker/reachability/SparseSccModelChecker.h" #include "src/storage/parameters.h" +#include "src/models/Dtmc.h" +#include "src/properties/prctl/PrctlFilter.h" + +//std::tuple>, boost::optional>, boost::optional, boost::optional, boost::optional> computeReachabilityProbability(storm::models::Dtmc const& dtmc, std::shared_ptr> const& filterFormula) { +// // The first thing we need to do is to make sure the formula is of the correct form and - if so - extract +// // the bitvector representation of the atomic propositions. +// +// std::shared_ptr> stateFormula = std::dynamic_pointer_cast>(filterFormula->getChild()); +// std::shared_ptr> pathFormula; +// boost::optional threshold; +// boost::optional strict; +// if (stateFormula != nullptr) { +// std::shared_ptr> probabilisticBoundFormula = std::dynamic_pointer_cast>(stateFormula); +// STORM_LOG_THROW(probabilisticBoundFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted."); +// STORM_LOG_THROW(probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS_EQUAL || probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties with upper probability bounds are permitted."); +// +// threshold = probabilisticBoundFormula->getBound(); +// strict = probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS; +// pathFormula = probabilisticBoundFormula->getChild(); +// } else { +// pathFormula = std::dynamic_pointer_cast>(filterFormula->getChild()); +// } +// +// STORM_LOG_THROW(pathFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted."); +// +// std::shared_ptr> untilFormula = std::dynamic_pointer_cast>(pathFormula); +// std::shared_ptr> phiStateFormula; +// std::shared_ptr> psiStateFormula; +// if (untilFormula != nullptr) { +// phiStateFormula = untilFormula->getLeft(); +// psiStateFormula = untilFormula->getRight(); +// } else { +// std::shared_ptr> eventuallyFormula = std::dynamic_pointer_cast>(pathFormula); +// STORM_LOG_THROW(eventuallyFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted."); +// phiStateFormula = std::shared_ptr>(new storm::properties::prctl::Ap("true")); +// psiStateFormula = eventuallyFormula->getChild(); +// } +// +// // Now we need to make sure the formulas defining the phi and psi states are just labels. +// std::shared_ptr> phiStateFormulaApFormula = std::dynamic_pointer_cast>(phiStateFormula); +// std::shared_ptr> psiStateFormulaApFormula = std::dynamic_pointer_cast>(psiStateFormula); +// STORM_LOG_THROW(phiStateFormulaApFormula.get() != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *phiStateFormula << " for parametric model checking. Note that only atomic propositions are admitted in that position."); +// STORM_LOG_THROW(psiStateFormulaApFormula.get() != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *psiStateFormula << " for parametric model checking. Note that only atomic propositions are admitted in that position."); +// +// // Now retrieve the appropriate bitvectors from the atomic propositions. +// storm::storage::BitVector phiStates = phiStateFormulaApFormula->getAp() != "true" ? dtmc.getLabeledStates(phiStateFormulaApFormula->getAp()) : storm::storage::BitVector(dtmc.getNumberOfStates(), true); +// storm::storage::BitVector psiStates = dtmc.getLabeledStates(psiStateFormulaApFormula->getAp()); +// +// // Do some sanity checks to establish some required properties. +// STORM_LOG_THROW(dtmc.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); +// +// // Then, compute the subset of states that has a probability of 0 or 1, respectively. +// std::pair statesWithProbability01 = storm::utility::graph::performProb01(dtmc, phiStates, psiStates); +// 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 (dtmc.getInitialStates().isDisjointFrom(maybeStates)) { +// STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step."); +// return statesWithProbability0.get(*dtmc.getInitialStates().begin()) ? storm::utility::constantZero() : storm::utility::constantOne(); +// } +// +// // 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(dtmc.getTransitionMatrix(), dtmc.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 = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1); +// +// // Determine the set of initial states of the sub-DTMC. +// storm::storage::BitVector newInitialStates = dtmc.getInitialStates() % maybeStates; +// +// // We then build the submatrix that only has the transitions of the maybe states. +// storm::storage::SparseMatrix submatrix = dtmc.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates); +// +// // To be able to apply heuristics later, we now determine the distance of each state to the initial state. +// std::vector> stateQueue; +// stateQueue.reserve(submatrix.getRowCount()); +// storm::storage::BitVector statesInQueue(submatrix.getRowCount()); +// std::vector distances(submatrix.getRowCount()); +// +// storm::storage::sparse::state_type currentPosition = 0; +// for (auto const& initialState : newInitialStates) { +// stateQueue.emplace_back(initialState, 0); +// statesInQueue.set(initialState); +// } +// +// // Perform a BFS. +// while (currentPosition < stateQueue.size()) { +// std::pair const& stateDistancePair = stateQueue[currentPosition]; +// distances[stateDistancePair.first] = stateDistancePair.second; +// +// for (auto const& successorEntry : submatrix.getRow(stateDistancePair.first)) { +// if (!statesInQueue.get(successorEntry.getColumn())) { +// stateQueue.emplace_back(successorEntry.getColumn(), stateDistancePair.second + 1); +// statesInQueue.set(successorEntry.getColumn()); +// } +// } +// ++currentPosition; +// } +// +// storm::modelchecker::reachability::SparseSccModelChecker modelchecker; +// +// return std::make_tuple(modelchecker.computeReachabilityProbability(submatrix, oneStepProbabilities, submatrix.transpose(), newInitialStates, phiStates, psiStates, distances),submatrix, oneStepProbabilities, newInitialStates, threshold, strict); +//} + /*! * Main entry point of the executable storm. */ int main(const int argc, const char** argv) { - try { +// try { storm::utility::cli::setUp(); storm::utility::cli::printHeader(argc, argv); bool optionsCorrect = storm::utility::cli::parseOptions(argc, argv); @@ -27,35 +138,57 @@ int main(const int argc, const char** argv) { // From this point on we are ready to carry out the actual computations. // Program Translation Time Measurement, Start std::chrono::high_resolution_clock::time_point programTranslationStart = std::chrono::high_resolution_clock::now(); - + // First, we build the model using the given symbolic model description and constant definitions. std::string const& programFile = storm::settings::generalSettings().getSymbolicModelFilename(); std::string const& constants = storm::settings::generalSettings().getConstantDefinitionString(); storm::prism::Program program = storm::parser::PrismParser::parse(programFile); std::shared_ptr> model = storm::adapters::ExplicitModelAdapter::translateProgram(program, constants); - - model->printModelInformationToStream(std::cout); - + // Program Translation Time Measurement, End std::chrono::high_resolution_clock::time_point programTranslationEnd = std::chrono::high_resolution_clock::now(); std::cout << "Parsing and translating the model took " << std::chrono::duration_cast(programTranslationEnd - programTranslationStart).count() << "ms." << std::endl << std::endl; - + std::shared_ptr> dtmc = model->as>(); + + // Perform bisimulation minimization if requested. + if (storm::settings::generalSettings().isBisimulationSet()) { + storm::storage::DeterministicModelStrongBisimulationDecomposition bisimulationDecomposition(*dtmc, true); + dtmc = bisimulationDecomposition.getQuotient()->as>(); + + dtmc->printModelInformationToStream(std::cout); + } + assert(dtmc); storm::modelchecker::reachability::CollectConstraints constraintCollector; constraintCollector(*dtmc); - - - storm::modelchecker::reachability::SparseSccModelChecker modelChecker; STORM_LOG_THROW(storm::settings::generalSettings().isPctlPropertySet(), storm::exceptions::InvalidSettingsException, "Unable to perform model checking without a property."); std::shared_ptr> filterFormula = storm::parser::PrctlParser::parsePrctlFormula(storm::settings::generalSettings().getPctlProperty()); - - storm::RationalFunction value = modelChecker.computeReachabilityProbability(*dtmc, filterFormula); - STORM_PRINT_AND_LOG(std::endl << "computed value " << value << std::endl); - + + storm::modelchecker::reachability::SparseSccModelChecker modelchecker; + + storm::RationalFunction valueFunction = modelchecker.computeReachabilityProbability(*dtmc, filterFormula); +// STORM_PRINT_AND_LOG(std::endl << "Result: (" << carl::computePolynomial(valueFunction.nominator()) << ") / (" << carl::computePolynomial(valueFunction.denominator()) << ")" << std::endl); + STORM_PRINT_AND_LOG(std::endl << "Result: (" << valueFunction.nominator() << ") / (" << valueFunction.denominator() << ")" << std::endl); + STORM_PRINT_AND_LOG(std::endl << "Result: " << valueFunction << std::endl); + +// // Perform bisimulation minimization if requested. +// if (storm::settings::generalSettings().isBisimulationSet()) { +// storm::storage::DeterministicModelStrongBisimulationDecomposition bisimulationDecomposition(*dtmc, true); +// dtmc = bisimulationDecomposition.getQuotient()->as>(); +// +// dtmc->printModelInformationToStream(std::cout); +// } + +// storm::RationalFunction valueFunction2 = modelchecker.computeReachabilityProbability(*dtmc, filterFormula); +// STORM_PRINT_AND_LOG(std::endl << "computed value2 " << valueFunction2 << std::endl); +// +// storm::RationalFunction diff = storm::utility::simplify(valueFunction - valueFunction2); +// STORM_PRINT_AND_LOG(std::endl << "difference: " << diff << std::endl); + // Get variables from parameter definitions in prism program. std::set parameters; for(auto constant : program.getConstants()) @@ -67,27 +200,26 @@ int main(const int argc, const char** argv) { parameters.insert(p); } } - // - STORM_LOG_ASSERT(parameters == value.gatherVariables(), "Parameters in result and program definition do not coincide."); + + STORM_LOG_ASSERT(parameters == valueFunction.gatherVariables(), "Parameters in result and program definition do not coincide."); if(storm::settings::parametricSettings().exportResultToFile()) { - storm::utility::exportParametricMcResult(value, constraintCollector); + storm::utility::exportParametricMcResult(valueFunction, constraintCollector); } - if(storm::settings::parametricSettings().exportToSmt2File()) { - storm::modelchecker::reachability::DirectEncoding dec; - //storm::utility::exportStreamToFile(dec.encodeAsSmt2(modelChecker.getMatrix(), parameters,)); - } - - +// if (storm::settings::parametricSettings().exportToSmt2File() && std::get<1>(result) && std::get<2>(result) && std::get<3>(result) && std::get<4>(result) && std::get<5>(result)) { +// storm::modelchecker::reachability::DirectEncoding dec; +// storm::utility::exportStringStreamToFile(dec.encodeAsSmt2(std::get<1>(result).get(), std::get<2>(result).get(), parameters, std::get<3>(result).get(), carl::rationalize(std::get<4>(result).get()), std::get<5>(result).get()), "out.smt"); +// } + // All operations have now been performed, so we clean up everything and terminate. storm::utility::cli::cleanUp(); return 0; - } catch (storm::exceptions::BaseException const& exception) { - STORM_LOG_ERROR("An exception caused StoRM to terminate. The message of the exception is: " << exception.what()); - } catch (std::exception const& exception) { - STORM_LOG_ERROR("An unexpected exception occurred and caused StoRM to terminate. The message of this exception is: " << exception.what()); - } +// } catch (storm::exceptions::BaseException const& exception) { +// STORM_LOG_ERROR("An exception caused StoRM to terminate. The message of the exception is: " << exception.what()); +// } catch (std::exception const& exception) { +// STORM_LOG_ERROR("An unexpected exception occurred and caused StoRM to terminate. The message of this exception is: " << exception.what()); +// } } //#include @@ -114,33 +246,33 @@ int main(const int argc, const char** argv) { // //std::string ParametricStormEntryPoint::reachabilityToSmt2(std::string const& label) //{ -// +// // storm::storage::BitVector phiStates(mModel->getNumberOfStates(), true); // storm::storage::BitVector initStates = mModel->getInitialStates(); // storm::storage::BitVector targetStates = mModel->getLabeledStates(label); -// +// // std::shared_ptr> dtmc = mModel->as>(); // // 1. make target states absorbing. // dtmc->makeAbsorbing(targetStates); // // 2. throw away anything which does not add to the reachability probability. // // 2a. remove non productive states // storm::storage::BitVector productiveStates = utility::graph::performProbGreater0(*dtmc, dtmc->getBackwardTransitions(), phiStates, targetStates); -// // 2b. calculate set of states wich +// // 2b. calculate set of states wich // storm::storage::BitVector almostSurelyReachingTargetStates = ~utility::graph::performProbGreater0(*dtmc, dtmc->getBackwardTransitions(), phiStates, ~productiveStates); // // 2c. Make such states also target states. // dtmc->makeAbsorbing(almostSurelyReachingTargetStates); -// // 2d. throw away non reachable states +// // 2d. throw away non reachable states // storm::storage::BitVector reachableStates = utility::graph::performProbGreater0(*dtmc, dtmc->getTransitionMatrix(), phiStates, initStates); // storm::storage::BitVector bv = productiveStates & reachableStates; // dtmc->getStateLabeling().addAtomicProposition("__targets__", targetStates | almostSurelyReachingTargetStates); // models::Dtmc subdtmc = dtmc->getSubDtmc(bv); -// +// // phiStates = storm::storage::BitVector(subdtmc.getNumberOfStates(), true); // initStates = subdtmc.getInitialStates(); // targetStates = subdtmc.getLabeledStates("__targets__"); // storm::storage::BitVector deadlockStates(phiStates); // deadlockStates.set(subdtmc.getNumberOfStates()-1,false); -// +// // // Search for states with only one non-deadlock successor. // std::map> chainedStates; // StateId nrStates = subdtmc.getNumberOfStates(); @@ -192,8 +324,8 @@ int main(const int argc, const char** argv) { // eliminatedStates.set(it->first, true); // } // } -// -// +// +// // for(auto chainedState : chainedStates) // { // if(!eliminatedStates[chainedState.first]) @@ -201,15 +333,15 @@ int main(const int argc, const char** argv) { // std::cout << chainedState.first << " -- " << chainedState.second.probability() << " --> " << chainedState.second.targetState() << std::endl; // } // } -// +// // storage::StronglyConnectedComponentDecomposition sccs(subdtmc); // std::cout << sccs << std::endl; // // modelchecker::reachability::DirectEncoding dec; // std::vector parameters; - + // return dec.encodeAsSmt2(subdtmc, parameters, subdtmc.getLabeledStates("init"), subdtmc.getLabeledStates("__targets__"), mpq_class(1,2)); -// +// //} // // @@ -224,8 +356,8 @@ int main(const int argc, const char** argv) { // fstream << entry.reachabilityToSmt2(s->getOptionByLongName("reachability").getArgument(0).getValueAsString()); // fstream.close(); // } -// -// +// +// //} // //} diff --git a/src/utility/ConstantsComparator.cpp b/src/utility/ConstantsComparator.cpp new file mode 100644 index 000000000..fd5fc9383 --- /dev/null +++ b/src/utility/ConstantsComparator.cpp @@ -0,0 +1,148 @@ +#include "src/utility/ConstantsComparator.h" + +#include "src/storage/sparse/StateType.h" + +namespace storm { + namespace utility { + + template + ValueType one() { + return ValueType(1); + } + + template + ValueType zero() { + return ValueType(0); + } + + template + ValueType infinity() { + return std::numeric_limits::infinity(); + } + + template + ValueType simplify(ValueType value) { + // In the general case, we don't to anything here, but merely return the value. If something else is + // supposed to happen here, the templated function can be specialized for this particular type. + return std::forward(value); + } + + template + bool ConstantsComparator::isOne(ValueType const& value) const { + return value == one(); + } + + template + bool ConstantsComparator::isZero(ValueType const& value) const { + return value == zero(); + } + + template + bool ConstantsComparator::isEqual(ValueType const& value1, ValueType const& value2) const { + return value1 == value2; + } + + ConstantsComparator::ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) { + // Intentionally left empty. + } + + ConstantsComparator::ConstantsComparator(double precision) : precision(precision) { + // Intentionally left empty. + } + + bool ConstantsComparator::isOne(double const& value) const { + return std::abs(value - one()) <= precision; + } + + bool ConstantsComparator::isZero(double const& value) const { + return std::abs(value) <= precision; + } + + bool ConstantsComparator::isEqual(double const& value1, double const& value2) const { + return std::abs(value1 - value2) <= precision; + } + + bool ConstantsComparator::isConstant(double const& value) const { + return true; + } + +#ifdef PARAMETRIC_SYSTEMS + template<> + RationalFunction& simplify(RationalFunction& value) { + value.simplify(); + return value; + } + + template<> + RationalFunction&& simplify(RationalFunction&& value) { + value.simplify(); + return std::move(value); + } + + bool ConstantsComparator::isOne(storm::RationalFunction const& value) const { + return value.isOne(); + } + + bool ConstantsComparator::isZero(storm::RationalFunction const& value) const { + return value.isZero(); + } + + bool ConstantsComparator::isEqual(storm::RationalFunction const& value1, storm::RationalFunction const& value2) const { + return value1 == value2; + } + + bool ConstantsComparator::isConstant(storm::RationalFunction const& value) const { + return value.isConstant(); + } + + bool ConstantsComparator::isOne(storm::Polynomial const& value) const { + return value.isOne(); + } + + bool ConstantsComparator::isZero(storm::Polynomial const& value) const { + return value.isZero(); + } + + bool ConstantsComparator::isEqual(storm::Polynomial const& value1, storm::Polynomial const& value2) const { + return value1 == value2; + } + + bool ConstantsComparator::isConstant(storm::Polynomial const& value) const { + return value.isConstant(); + } + +#endif + + template + storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry) { + simplify(matrixEntry.getValue()); + return matrixEntry; + } + + template class ConstantsComparator; + template class ConstantsComparator; + template class ConstantsComparator; + + template double one(); + template double zero(); + template double infinity(); + + template RationalFunction one(); + template RationalFunction zero(); + + template Polynomial one(); + template Polynomial zero(); + + template double simplify(double value); + template RationalFunction simplify(RationalFunction value); + + template double& simplify(double& value); + template RationalFunction& simplify(RationalFunction& value); + + template double&& simplify(double&& value); + template RationalFunction&& simplify(RationalFunction&& value); + + template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); + template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); + } +} \ No newline at end of file diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h new file mode 100644 index 000000000..37e27ea46 --- /dev/null +++ b/src/utility/ConstantsComparator.h @@ -0,0 +1,104 @@ +#ifndef STORM_UTILITY_CONSTANTSCOMPARATOR_H_ +#define STORM_UTILITY_CONSTANTSCOMPARATOR_H_ + +#ifdef max +# undef max +#endif + +#ifdef min +# undef min +#endif + +#include +#include + +#include "src/settings/SettingsManager.h" +#include "src/storage/SparseMatrix.h" + +namespace storm { + namespace utility { + + template + ValueType one(); + + template + ValueType zero(); + + template + ValueType infinity(); + + template + ValueType simplify(ValueType value); + + // A class that can be used for comparing constants. + template + class ConstantsComparator { + public: + bool isOne(ValueType const& value) const; + + bool isZero(ValueType const& value) const; + + bool isEqual(ValueType const& value1, ValueType const& value2) const; + }; + + // For doubles we specialize this class and consider the comparison modulo some predefined precision. + template<> + class ConstantsComparator { + public: + ConstantsComparator(); + + ConstantsComparator(double precision); + + bool isOne(double const& value) const; + + bool isZero(double const& value) const; + + bool isEqual(double const& value1, double const& value2) const; + + bool isConstant(double const& value) const; + + private: + // The precision used for comparisons. + double precision; + }; + +#ifdef PARAMETRIC_SYSTEMS + template<> + RationalFunction& simplify(RationalFunction& value); + + template<> + RationalFunction&& simplify(RationalFunction&& value); + + template<> + class ConstantsComparator { + public: + bool isOne(storm::RationalFunction const& value) const; + + bool isZero(storm::RationalFunction const& value) const; + + bool isEqual(storm::RationalFunction const& value1, storm::RationalFunction const& value2) const; + + bool isConstant(storm::RationalFunction const& value) const; + }; + + template<> + class ConstantsComparator { + public: + bool isOne(storm::Polynomial const& value) const; + + bool isZero(storm::Polynomial const& value) const; + + bool isEqual(storm::Polynomial const& value1, storm::Polynomial const& value2) const; + + bool isConstant(storm::Polynomial const& value) const; + }; + + +#endif + + template + storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); + } +} + +#endif /* STORM_UTILITY_CONSTANTSCOMPARATOR_H_ */ \ No newline at end of file diff --git a/src/utility/cli.h b/src/utility/cli.h index 54da82cb4..227e09dfb 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -46,6 +46,10 @@ log4cplus::Logger printer; // Headers of adapters. #include "src/adapters/ExplicitModelAdapter.h" +// Headers for model processing. +#include "src/storage/NaiveDeterministicModelBisimulationDecomposition.h" +#include "src/storage/DeterministicModelStrongBisimulationDecomposition.h" + // Headers for counterexample generation. #include "src/counterexamples/MILPMinimalLabelSetGenerator.h" #include "src/counterexamples/SMTMinimalCommandSetGenerator.h" @@ -259,6 +263,22 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "No input model."); } + // Print some information about the model. + result->printModelInformationToStream(std::cout); + + if (settings.isBisimulationSet()) { + STORM_LOG_THROW(result->getType() == storm::models::DTMC, storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only compatible with DTMCs."); + std::shared_ptr> dtmc = result->template as>(); + + STORM_PRINT(std::endl << "Performing bisimulation minimization..." << std::endl); + storm::storage::DeterministicModelStrongBisimulationDecomposition bisimulationDecomposition(*dtmc, true); + + result = bisimulationDecomposition.getQuotient(); + + STORM_PRINT_AND_LOG(std::endl << "Model after minimization:" << std::endl); + result->printModelInformationToStream(std::cout); + } + return result; } @@ -269,7 +289,7 @@ namespace storm { std::shared_ptr> mdp = model->as>(); - // FIXME: re-parse the program. + // FIXME: do not re-parse the program. std::string const programFile = storm::settings::generalSettings().getSymbolicModelFilename(); std::string const constants = storm::settings::generalSettings().getConstantDefinitionString(); storm::prism::Program program = storm::parser::PrismParser::parse(programFile); @@ -298,15 +318,13 @@ namespace storm { // Start by parsing/building the model. std::shared_ptr> model = buildModel(); - // Print some information about the model. - model->printModelInformationToStream(std::cout); - // If we were requested to generate a counterexample, we now do so. if (settings.isCounterexampleSet()) { STORM_LOG_THROW(settings.isPctlPropertySet(), storm::exceptions::InvalidSettingsException, "Unable to generate counterexample without a property."); std::shared_ptr> filterFormula = storm::parser::PrctlParser::parsePrctlFormula(settings.getPctlProperty()); generateCounterexample(model, filterFormula->getChild()); } + } } diff --git a/src/utility/constants.h b/src/utility/constants.h index 699ee7368..60e601aaa 100644 --- a/src/utility/constants.h +++ b/src/utility/constants.h @@ -254,8 +254,6 @@ inline bool isZero(storm::Polynomial const& p) } #endif - - template inline bool isOne(T const& sum) { diff --git a/src/utility/export.h b/src/utility/export.h index 7d66bd0b3..35386e6a8 100644 --- a/src/utility/export.h +++ b/src/utility/export.h @@ -25,7 +25,7 @@ void exportParametricMcResult(const ValueType& mcresult, storm::modelchecker::re // todo add checks. filestream << "!Parameters: "; std::set vars = mcresult.gatherVariables(); - std::copy(vars.begin(), vars.end(), std::ostream_iterator(filestream, " ")); + std::copy(vars.begin(), vars.end(), std::ostream_iterator(filestream, ", ")); filestream << std::endl; filestream << "!Result: " << mcresult << std::endl; filestream << "!Well-formed Constraints: " << std::endl; @@ -35,11 +35,11 @@ void exportParametricMcResult(const ValueType& mcresult, storm::modelchecker::re filestream.close(); } -void exportStringStreamToFile(std::stringstream& stream, std::string filepath) { +void exportStringStreamToFile(std::string const& str, std::string filepath) { // todo add checks. std::ofstream filestream; filestream.open(filepath); - filestream << stream.str(); + filestream << str; filestream.close(); } diff --git a/src/utility/graph.h b/src/utility/graph.h index 33f637877..cbfad8758 100644 --- a/src/utility/graph.h +++ b/src/utility/graph.h @@ -204,6 +204,25 @@ namespace storm { return result; } + /*! + * Computes the sets of states that have probability 0 or 1, respectively, of satisfying phi until psi in a + * deterministic model. + * + * @param backwardTransitions The backward transitions of the model whose graph structure to search. + * @param phiStates The set of all states satisfying phi. + * @param psiStates The set of all states satisfying psi. + * @return A pair of bit vectors such that the first bit vector stores the indices of all states + * with probability 0 and the second stores all indices of states with probability 1. + */ + template + static std::pair performProb01(storm::storage::SparseMatrix backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) { + std::pair result; + result.first = performProbGreater0(backwardTransitions, phiStates, psiStates); + result.second = performProb1(backwardTransitions, phiStates, psiStates, result.first); + result.first.complement(); + return result; + } + /*! * Computes the sets of states that have probability greater 0 of satisfying phi until psi under at least * one possible resolution of non-determinism in a non-deterministic model. Stated differently,