From 7156a63b0fa818a4427987b0726cbed21f514c3e Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 6 Nov 2015 18:58:26 +0100 Subject: [PATCH] tried different approach for bisim for MDPs Former-commit-id: 92d56a462076debaf95d974c2f65e933d9652361 --- src/builder/DdPrismModelBuilder.cpp | 2 +- src/builder/ExplicitPrismModelBuilder.cpp | 2 +- src/storage/Distribution.cpp | 67 +++- src/storage/Distribution.h | 43 ++- .../BisimulationDecomposition.cpp | 12 + .../bisimulation/BisimulationDecomposition.h | 5 + ...ministicModelBisimulationDecomposition.cpp | 2 +- ...ministicModelBisimulationDecomposition.cpp | 296 +++++++++--------- ...erministicModelBisimulationDecomposition.h | 45 ++- src/utility/ConstantsComparator.cpp | 15 +- src/utility/ConstantsComparator.h | 6 + src/utility/storm.h | 30 +- 12 files changed, 330 insertions(+), 195 deletions(-) diff --git a/src/builder/DdPrismModelBuilder.cpp b/src/builder/DdPrismModelBuilder.cpp index 4f6ab3871..0a20bf38b 100644 --- a/src/builder/DdPrismModelBuilder.cpp +++ b/src/builder/DdPrismModelBuilder.cpp @@ -1020,7 +1020,7 @@ namespace storm { preparedProgram = preparedProgram.substituteConstants(); - STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl); + STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl); // Start by initializing the structure used for storing all information needed during the model generation. // In particular, this creates the meta variables used to encode the model. diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp index 5b62364b4..4f6117c9c 100644 --- a/src/builder/ExplicitPrismModelBuilder.cpp +++ b/src/builder/ExplicitPrismModelBuilder.cpp @@ -295,7 +295,7 @@ namespace storm { // Now that the program is fixed, we we need to substitute all constants with their concrete value. preparedProgram = preparedProgram.substituteConstants(); - STORM_LOG_DEBUG("Building representation of program :" << std::endl << preparedProgram << std::endl); + STORM_LOG_DEBUG("Building representation of program:" << std::endl << preparedProgram << std::endl); // Select the appropriate reward models (after the constants have been substituted). std::vector> selectedRewardModels; diff --git a/src/storage/Distribution.cpp b/src/storage/Distribution.cpp index ce42b6ab4..5dd19e50c 100644 --- a/src/storage/Distribution.cpp +++ b/src/storage/Distribution.cpp @@ -3,8 +3,14 @@ #include #include +#include "src/utility/macros.h" +#include "src/utility/constants.h" +#include "src/utility/ConstantsComparator.h" + #include "src/settings/SettingsManager.h" +#include "src/adapters/CarlAdapter.h" + namespace storm { namespace storage { @@ -14,7 +20,7 @@ namespace storm { } template - bool Distribution::operator==(Distribution const& other) const { + bool Distribution::equals(Distribution const& other, storm::utility::ConstantsComparator const& comparator) 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; @@ -28,7 +34,7 @@ namespace storm { if (first1->first != first2->first) { return false; } - if (std::abs(first1->second - first2->second) > 1e-6) { + if (!comparator.isEqual(first1->second, first2->second)) { return false; } } @@ -38,10 +44,29 @@ namespace storm { template void Distribution::addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability) { - if (this->distribution.find(state) == this->distribution.end()) { + auto it = this->distribution.find(state); + if (it == this->distribution.end()) { this->hash += static_cast(state); + this->distribution.emplace_hint(it, state, probability); + } else { + it->second += probability; + } + } + + template + void Distribution::removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator const& comparator) { + auto it = this->distribution.find(state); + STORM_LOG_ASSERT(it != this->distribution.end(), "Cannot remove probability, because the state is not in the support of the distribution."); + it->second -= probability; + if (comparator.isZero(it->second)) { + this->distribution.erase(it); } - this->distribution[state] += probability; + } + + template + void Distribution::shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator const& comparator) { + removeProbability(fromState, probability, comparator); + addProbability(toState, probability); } template @@ -68,7 +93,7 @@ namespace storm { 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; + ValueType scaleValue = storm::utility::one() / probabilityIterator->second; this->distribution.erase(probabilityIterator); for (auto& entry : this->distribution) { @@ -82,6 +107,11 @@ namespace storm { return this->hash ^ (this->distribution.size() << 8); } + template + std::size_t Distribution::size() const { + return this->distribution.size(); + } + template std::ostream& operator<<(std::ostream& out, Distribution const& distribution) { out << "{"; @@ -93,7 +123,34 @@ namespace storm { return out; } + template + bool Distribution::less(Distribution const& other, storm::utility::ConstantsComparator const& comparator) const { + if (this->size() != other.size()) { + return this->size() < other.size(); + } + auto firstIt = this->begin(); + auto firstIte = this->end(); + auto secondIt = other.begin(); + for (; firstIt != firstIte; ++firstIt, ++secondIt) { + // If the two blocks already differ, we can decide which distribution is smaller. + if (firstIt->first != secondIt->first) { + return firstIt->first < secondIt->first; + } + + // If the blocks are the same, but the probability differs, we can also decide which distribution is smaller. + if (!comparator.isEqual(firstIt->second, secondIt->second)) { + return comparator.isLess(firstIt->second, secondIt->second); + } + } + return false; + } + template class Distribution; template std::ostream& operator<<(std::ostream& out, Distribution const& distribution); + +#ifdef STORM_HAVE_CARL + template class Distribution; + template std::ostream& operator<<(std::ostream& out, Distribution const& distribution); +#endif } } diff --git a/src/storage/Distribution.h b/src/storage/Distribution.h index 010cf43bc..93330c009 100644 --- a/src/storage/Distribution.h +++ b/src/storage/Distribution.h @@ -8,6 +8,11 @@ #include "src/storage/sparse/StateType.h" namespace storm { + namespace utility { + template + class ConstantsComparator; + } + namespace storage { template @@ -28,7 +33,7 @@ namespace storm { * @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; + bool equals(Distribution const& other, storm::utility::ConstantsComparator const& comparator = storm::utility::ConstantsComparator()) const; /*! * Assigns the given state the given probability under this distribution. @@ -38,13 +43,34 @@ namespace storm { */ void addProbability(storm::storage::sparse::state_type const& state, ValueType const& probability); + /*! + * Removes the given probability mass of going to the given state. + * + * @param state The state for which to remove the probability. + * @param probability The probability to remove. + * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the + * entry is removed. + */ + void removeProbability(storm::storage::sparse::state_type const& state, ValueType const& probability, storm::utility::ConstantsComparator const& comparator = storm::utility::ConstantsComparator()); + + /*! + * Removes the probability mass from one state and adds it to another. + * + * @param fromState The state from which to take the probability mass. + * @param toState The state from which to which to add the probability mass. + * @param probability The probability mass to shift. + * @param comparator A comparator that is used to determine if the remaining probability is zero. If so, the + * entry is removed. + */ + void shiftProbability(storm::storage::sparse::state_type const& fromState, storm::storage::sparse::state_type const& toState, ValueType const& probability, storm::utility::ConstantsComparator const& comparator = storm::utility::ConstantsComparator()); + /*! * 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. * @@ -58,7 +84,7 @@ namespace storm { * @return The iterator past the elements in this distribution. */ iterator end(); - + /*! * Retrieves an iterator past the elements in this distribution. * @@ -82,16 +108,23 @@ namespace storm { */ std::size_t getHash() const; + /*! + * Retrieves the size of the distribution, i.e. the size of the support set. + */ + std::size_t size() const; + + bool less(Distribution const& other, storm::utility::ConstantsComparator const& comparator) 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 + // A hash value that is maintained to allow for quicker equality comparison between distributions. std::size_t hash; }; template - std::ostream& operator<<(std::ostream& out, Distribution const& distribution); + std::ostream& operator<<(std::ostream& out, Distribution const& distribution); } } diff --git a/src/storage/bisimulation/BisimulationDecomposition.cpp b/src/storage/bisimulation/BisimulationDecomposition.cpp index 86ac56378..0dc8369fa 100644 --- a/src/storage/bisimulation/BisimulationDecomposition.cpp +++ b/src/storage/bisimulation/BisimulationDecomposition.cpp @@ -4,6 +4,7 @@ #include "src/models/sparse/Dtmc.h" #include "src/models/sparse/Ctmc.h" +#include "src/models/sparse/Mdp.h" #include "src/models/sparse/StandardRewardModel.h" #include "src/storage/bisimulation/DeterministicBlockData.h" @@ -167,8 +168,12 @@ namespace storm { } else { this->initializeLabelBasedPartition(); } + std::cout << "initial partition is " << std::endl; + this->partition.print(); std::chrono::high_resolution_clock::duration initialPartitionTime = std::chrono::high_resolution_clock::now() - initialPartitionStart; + this->initialize(); + std::chrono::high_resolution_clock::time_point refinementStart = std::chrono::high_resolution_clock::now(); this->performPartitionRefinement(); std::chrono::high_resolution_clock::duration refinementTime = std::chrono::high_resolution_clock::now() - refinementStart; @@ -274,6 +279,11 @@ namespace storm { } } + template + void BisimulationDecomposition::initialize() { + // Intentionally left empty. + } + template void BisimulationDecomposition::extractDecompositionBlocks() { // Now move the states from the internal partition into their final place in the decomposition. We do so in @@ -290,10 +300,12 @@ namespace storm { template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; + template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; #ifdef STORM_HAVE_CARL template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; + template class BisimulationDecomposition, bisimulation::DeterministicBlockData>; #endif } } diff --git a/src/storage/bisimulation/BisimulationDecomposition.h b/src/storage/bisimulation/BisimulationDecomposition.h index 280737797..a7bd575b2 100644 --- a/src/storage/bisimulation/BisimulationDecomposition.h +++ b/src/storage/bisimulation/BisimulationDecomposition.h @@ -185,6 +185,11 @@ namespace storm { */ virtual void initializeMeasureDrivenPartition(); + /*! + * A function that can initialize auxiliary data structures. It is called after initializing the initial partition. + */ + virtual void initialize(); + /*! * Computes the set of states with probability 0/1 for satisfying phi until psi. This is used for the measure * driven initial partition. diff --git a/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp index e3867f2de..c9672cd47 100644 --- a/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp +++ b/src/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp @@ -172,7 +172,7 @@ namespace storm { split |= this->partition.splitBlock(*blockToRefineProbabilistically, [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { - return getProbabilityToSplitter(state1) < getProbabilityToSplitter(state2); + return this->comparator.isLess(getProbabilityToSplitter(state1), getProbabilityToSplitter(state2)); }, [&splitterQueue] (Block& block) { splitterQueue.emplace_back(&block); block.data().setSplitter(); diff --git a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp index e519e1b2f..1b16cd6b4 100644 --- a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp +++ b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp @@ -5,19 +5,21 @@ #include "src/utility/graph.h" +#include "src/utility/macros.h" #include "src/exceptions/IllegalFunctionCallException.h" +#include "src/adapters/CarlAdapter.h" + namespace storm { namespace storage { - + using namespace bisimulation; template - NondeterministicModelBisimulationDecomposition::NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition::Options const& options) : BisimulationDecomposition(model, model.getTransitionMatrix().transpose(false), options), choiceToStateMapping(model.getNumberOfChoices()), probabilitiesToCurrentSplitter(model.getNumberOfChoices(), storm::utility::zero()) { + NondeterministicModelBisimulationDecomposition::NondeterministicModelBisimulationDecomposition(ModelType const& model, typename BisimulationDecomposition::Options const& options) : BisimulationDecomposition(model, model.getTransitionMatrix().transpose(false), options), choiceToStateMapping(model.getNumberOfChoices()), quotientDistributions(model.getNumberOfChoices()), orderedQuotientDistributions(model.getNumberOfChoices()) { STORM_LOG_THROW(options.type == BisimulationType::Strong, storm::exceptions::IllegalFunctionCallException, "Weak bisimulation is currently not supported for nondeterministic models."); - this->createChoiceToStateMapping(); } - + template std::pair NondeterministicModelBisimulationDecomposition::getStatesWithProbability01() { STORM_LOG_THROW(static_cast(this->options.optimalityType), storm::exceptions::IllegalFunctionCallException, "Can only compute states with probability 0/1 with an optimization direction (min/max)."); @@ -27,7 +29,13 @@ namespace storm { return storm::utility::graph::performProb01Max(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(), this->model.getBackwardTransitions(), this->options.phiStates.get(), this->options.psiStates.get()); } } - + + template + void NondeterministicModelBisimulationDecomposition::initialize() { + this->createChoiceToStateMapping(); + this->initializeQuotientDistributions(); + } + template void NondeterministicModelBisimulationDecomposition::createChoiceToStateMapping() { std::vector nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices(); @@ -39,192 +47,172 @@ namespace storm { } template - void NondeterministicModelBisimulationDecomposition::buildQuotient() { - STORM_LOG_ASSERT(false, "Not yet implemented"); - } - - template - bool NondeterministicModelBisimulationDecomposition::possiblyNeedsRefinement(bisimulation::Block const& predecessorBlock) const { - return predecessorBlock.getNumberOfStates() > 1 && !predecessorBlock.data().absorbing(); - } - - template - void NondeterministicModelBisimulationDecomposition::clearProbabilitiesToSplitter(storm::storage::sparse::state_type state) { + void NondeterministicModelBisimulationDecomposition::initializeQuotientDistributions() { std::vector nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices(); - for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) { - probabilitiesToCurrentSplitter[choice] = storm::utility::zero(); + for (auto choice = 0; choice < nondeterministicChoiceIndices.back(); ++choice) { + for (auto entry : this->model.getTransitionMatrix().getRow(choice)) { + if (!this->comparator.isZero(entry.getValue())) { + this->quotientDistributions[choice].addProbability(this->partition.getBlock(entry.getColumn()).getId(), entry.getValue()); + } + } + orderedQuotientDistributions[choice] = &this->quotientDistributions[choice]; } - } - - template - void NondeterministicModelBisimulationDecomposition::increaseProbabilityToSplitter(storm::storage::sparse::state_type state, uint_fast64_t choice, bisimulation::Block const& predecessorBlock, ValueType const& value) { - storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(state); - // If the position of the state is to the right of marker1, we have not seen it before. That means we need to - // clear the probability associated with all choices of the given state. - if (predecessorPosition >= predecessorBlock.data().marker1()) { - clearProbabilitiesToSplitter(state); + for (auto state = 0; state < this->model.getNumberOfStates(); ++state) { + updateOrderedQuotientDistributions(state); } - - // Now increase the probability of the given choice of the given state. - probabilitiesToCurrentSplitter[choice] += value; - } - - template - void NondeterministicModelBisimulationDecomposition::moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock) { - this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1())); - predecessorBlock.data().incrementMarker1(); } - template - void NondeterministicModelBisimulationDecomposition::moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock) { - this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker2())); - predecessorBlock.data().incrementMarker2(); + template + void NondeterministicModelBisimulationDecomposition::updateOrderedQuotientDistributions(storm::storage::sparse::state_type state) { + std::vector nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices(); + std::sort(this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state], this->orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state + 1], + [this] (storm::storage::Distribution const* dist1, storm::storage::Distribution const* dist2) { + return dist1->less(*dist2, this->comparator); + }); } - template - void NondeterministicModelBisimulationDecomposition::moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip) { - storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor); - - // If the predecessors of the given predecessor were already explored, we can move it easily. - if (predecessorPosition <= currentPositionInSplitter + elementsToSkip) { - this->partition.swapStates(predecessor, this->partition.getState(predecessorBlock.data().marker1())); - predecessorBlock.data().incrementMarker1(); - } else { - // Otherwise, we need to move the predecessor, but we need to make sure that we explore its - // predecessors later. We do this by moving it to a range at the beginning of the block that will hold - // all predecessors in the splitter whose predecessors have yet to be explored. - if (predecessorBlock.data().marker2() == predecessorBlock.data().marker1()) { - this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition); - this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1); - } else { - this->partition.swapStatesAtPositions(predecessorBlock.data().marker2(), predecessorPosition); - this->partition.swapStatesAtPositions(predecessorPosition, predecessorBlock.data().marker1()); - this->partition.swapStatesAtPositions(predecessorPosition, currentPositionInSplitter + elementsToSkip + 1); - } - - // Since we had to move an already explored state to the right of the current position, - ++elementsToSkip; - predecessorBlock.data().incrementMarker1(); - predecessorBlock.data().incrementMarker2(); - } + template + void NondeterministicModelBisimulationDecomposition::buildQuotient() { + std::cout << "Found " << this->partition.size() << " blocks" << std::endl; + this->partition.print(); + STORM_LOG_ASSERT(false, "Not yet implemented"); } template - void NondeterministicModelBisimulationDecomposition::insertIntoPredecessorList(bisimulation::Block& predecessorBlock, std::list*>& predecessorBlocks) { - // Insert the block into the list of blocks to refine (if that has not already happened). - if (!predecessorBlock.data().needsRefinement()) { - predecessorBlocks.emplace_back(&predecessorBlock); - predecessorBlock.data().setNeedsRefinement(); - } + bool NondeterministicModelBisimulationDecomposition::possiblyNeedsRefinement(bisimulation::Block const& block) const { + return block.getNumberOfStates() > 1 && !block.data().absorbing(); } - template - void NondeterministicModelBisimulationDecomposition::exploreRemainingStatesOfSplitter(bisimulation::Block& splitter, std::list*>& predecessorBlocks) { - for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.begin(splitter) + (splitter.data().marker2() - splitter.getBeginIndex()); splitterIt != splitterIte; ++splitterIt) { - storm::storage::sparse::state_type currentState = *splitterIt; - - for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) { - storm::storage::sparse::state_type choice = predecessorEntry.getColumn(); - storm::storage::sparse::state_type predecessor = choiceToStateMapping[choice]; - Block& predecessorBlock = this->partition.getBlock(predecessor); - - // If the block does not need to be refined, we skip it. - if (!possiblyNeedsRefinement(predecessorBlock)) { + template + void NondeterministicModelBisimulationDecomposition::updateQuotientDistributionsOfPredecessors(Block const& newBlock, Block const& oldBlock, std::deque*>& splitterQueue) { + uint_fast64_t lastState = 0; + bool lastStateInitialized = false; + + for (auto stateIt = this->partition.begin(newBlock), stateIte = this->partition.end(newBlock); stateIt != stateIte; ++stateIt) { + for (auto predecessorEntry : this->backwardTransitions.getRow(*stateIt)) { + if (this->comparator.isZero(predecessorEntry.getValue())) { continue; } - // We keep track of the probability of the predecessor moving to the splitter. - increaseProbabilityToSplitter(predecessor, choice, predecessorBlock, predecessorEntry.getValue()); + storm::storage::sparse::state_type predecessorChoice = predecessorEntry.getColumn(); + storm::storage::sparse::state_type predecessorState = choiceToStateMapping[predecessorChoice]; + Block& predecessorBlock = this->partition.getBlock(predecessorState); + + // If the predecessor block is not marked as to-refined, we do so now. + if (!predecessorBlock.data().splitter()) { + predecessorBlock.data().setSplitter(); + splitterQueue.push_back(&predecessorBlock); + } - // Only move the state if it has not been seen as a predecessor before. - storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor); - if (predecessorPosition >= predecessorBlock.data().marker1()) { - moveStateToMarker1(predecessor, predecessorBlock); + if (lastStateInitialized) { + // If we have skipped to the choices of the next state, we need to repair the order of the + // distributions for the last state. + if (lastState != predecessorState) { + updateOrderedQuotientDistributions(lastState); + lastState = predecessorState; + } + } else { + lastStateInitialized = true; + lastState = choiceToStateMapping[predecessorChoice]; } - insertIntoPredecessorList(predecessorBlock, predecessorBlocks); + // Now shift the probability from this transition from the old block to the new one. + this->quotientDistributions[predecessorChoice].shiftProbability(oldBlock.getId(), newBlock.getId(), predecessorEntry.getValue()); } } - // Finally, we can reset the second marker. - splitter.data().setMarker2(splitter.getBeginIndex()); + if (lastStateInitialized) { + updateOrderedQuotientDistributions(lastState); + } } template - void NondeterministicModelBisimulationDecomposition::refinePredecessorBlocksOfSplitter(std::list*> const& predecessorBlocks, std::deque*>& splitterQueue) { - // TODO + bool NondeterministicModelBisimulationDecomposition::checkQuotientDistributions() const { + } template - void NondeterministicModelBisimulationDecomposition::refinePartitionBasedOnSplitter(bisimulation::Block& splitter, std::deque*>& splitterQueue) { - // The outline of the refinement is as follows. - // - // We iterate over all states of the splitter and determine for each predecessor the state the probability - // entering the splitter. These probabilities are written to a member vector so that after the iteration - // process we have the probabilities of all predecessors of the splitter of entering the splitter in one - // step. To directly separate the states having a transition into the splitter from the ones who do not, - // we move the states to certain locations. That is, on encountering a predecessor of the splitter, it is - // moved to the beginning of its block. If the predecessor is in the splitter itself, we have to be a bit - // careful about where to move states. - // - // After this iteration, there may be states of the splitter whose predecessors have not yet been explored, - // so this needs to be done now. - // - // Finally, we use the information obtained in the first part for the actual splitting process in which all - // predecessor blocks of the splitter are split based on the probabilities computed earlier. - std::list*> predecessorBlocks; - storm::storage::sparse::state_type currentPosition = splitter.getBeginIndex(); - bool splitterIsPredecessorBlock = false; - for (auto splitterIt = this->partition.begin(splitter), splitterIte = this->partition.end(splitter); splitterIt != splitterIte; ++splitterIt, ++currentPosition) { - storm::storage::sparse::state_type currentState = *splitterIt; - - uint_fast64_t elementsToSkip = 0; - for (auto const& predecessorEntry : this->backwardTransitions.getRow(currentState)) { - storm::storage::sparse::state_type choice = predecessorEntry.getColumn(); - storm::storage::sparse::state_type predecessor = choiceToStateMapping[choice]; - storm::storage::sparse::state_type predecessorPosition = this->partition.getPosition(predecessor); - Block& predecessorBlock = this->partition.getBlock(predecessor); - - // If the block does not need to be refined, we skip it. - if (!possiblyNeedsRefinement(predecessorBlock)) { - continue; - } - - // We keep track of the probability of the predecessor moving to the splitter. - increaseProbabilityToSplitter(predecessor, choice, predecessorBlock, predecessorEntry.getValue()); - - // We only need to move the predecessor if its not already known as a predecessor already. - if (predecessorPosition >= predecessorBlock.data().marker1()) { - // If the predecessor block is not the splitter, we can move the state easily. - if (predecessorBlock != splitter) { - moveStateToMarker1(predecessor, predecessorBlock); - } else { - // If the predecessor is in the splitter, we need to be a bit more careful. - splitterIsPredecessorBlock = true; - moveStateInSplitter(predecessor, predecessorBlock, currentPosition, elementsToSkip); - } - - insertIntoPredecessorList(predecessorBlock, predecessorBlocks); - } + bool NondeterministicModelBisimulationDecomposition::splitBlockAccordingToCurrentQuotientDistributions(Block& block, std::deque*>& splitterQueue) { + bool split = this->partition.splitBlock(block, + [this] (storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { + return quotientDistributionsLess(state1, state2); + }, + [this, &block, &splitterQueue] (Block& newBlock) { + updateQuotientDistributionsOfPredecessors(newBlock, block, splitterQueue); + }); + + // The quotient distributions of the predecessors of block do not need to be updated, since the probability + // will go to the block with the same id as before. + +// std::cout << "partition after split: " << std::endl; +// this->partition.print(); + + this->checkQuotientDistributions(); + + return split; + } + + template + bool NondeterministicModelBisimulationDecomposition::quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2) { + STORM_LOG_TRACE("Comparing the quotient distributions of state " << state1 << " and " << state2 << "."); + std::vector nondeterministicChoiceIndices = this->model.getTransitionMatrix().getRowGroupIndices(); + +// std::cout << "state " << state1 << std::endl; +// for (int c = nondeterministicChoiceIndices[state1]; c < nondeterministicChoiceIndices[state1 + 1]; ++c) { +// std::cout << quotientDistributions[c] << std::endl; +// } +// +// std::cout << "state " << state2 << std::endl; +// for (int c = nondeterministicChoiceIndices[state2]; c < nondeterministicChoiceIndices[state2 + 1]; ++c) { +// std::cout << quotientDistributions[c] << std::endl; +// } + + auto firstIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1]; + auto firstIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state1 + 1]; + auto secondIt = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2]; + auto secondIte = orderedQuotientDistributions.begin() + nondeterministicChoiceIndices[state2 + 1]; + + for (; firstIt != firstIte && secondIt != secondIte; ++firstIt, ++secondIt) { + // If the current distributions are in a less-than relationship, we can return a result. + if ((*firstIt)->less(**secondIt, this->comparator)) { + return true; + } else if ((*secondIt)->less(**firstIt, this->comparator)) { + return false; } - // If, as a consequence of shifting states, we need to skip some elements, do so now. - splitterIt += elementsToSkip; - currentPosition += elementsToSkip; + // If the distributions matched, we need to advance both distribution iterators to the next distribution + // that is larger. + while (firstIt != firstIte && std::next(firstIt) != firstIte && !(*firstIt)->less(**std::next(firstIt), this->comparator)) { + ++firstIt; + } + while (secondIt != secondIte && std::next(secondIt) != secondIte && !(*secondIt)->less(**std::next(secondIt), this->comparator)) { + ++secondIt; + } } - // If the splitter was a predecessor block of itself, we potentially need to explore some states that have - // not been explored yet. - if (splitterIsPredecessorBlock) { - exploreRemainingStatesOfSplitter(splitter, predecessorBlocks); + if (firstIt == firstIte && secondIt != secondIte) { + return true; + } + return false; + } + + template + void NondeterministicModelBisimulationDecomposition::refinePartitionBasedOnSplitter(bisimulation::Block& splitter, std::deque*>& splitterQueue) { + if (!possiblyNeedsRefinement(splitter)) { + return; } - // Finally, we split the block based on the precomputed probabilities and the chosen bisimulation type. - refinePredecessorBlocksOfSplitter(predecessorBlocks, splitterQueue); + STORM_LOG_TRACE("Refining block " << splitter.getId()); + + splitBlockAccordingToCurrentQuotientDistributions(splitter, splitterQueue); } template class NondeterministicModelBisimulationDecomposition>; +#ifdef STORM_HAVE_CARL + template class NondeterministicModelBisimulationDecomposition>; +#endif + } } diff --git a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h index eb315dce6..827dc39cc 100644 --- a/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h +++ b/src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h @@ -4,6 +4,8 @@ #include "src/storage/bisimulation/BisimulationDecomposition.h" #include "src/storage/bisimulation/DeterministicBlockData.h" +#include "src/storage/Distribution.h" + namespace storm { namespace utility { template class ConstantsComparator; @@ -37,44 +39,41 @@ namespace storm { virtual void refinePartitionBasedOnSplitter(bisimulation::Block& splitter, std::deque*>& splitterQueue) override; + virtual void initialize() override; + private: // Creates the mapping from the choice indices to the states. void createChoiceToStateMapping(); - // Retrieves whether the given predecessor of the splitters possibly needs refinement. - bool possiblyNeedsRefinement(bisimulation::Block const& predecessorBlock) const; - - // Increases the probability of moving to the current splitter for the given choice of the given state. - void increaseProbabilityToSplitter(storm::storage::sparse::state_type state, uint_fast64_t choice, bisimulation::Block const& predecessorBlock, ValueType const& value); - - // Clears the probabilities of all choices of the given state. - void clearProbabilitiesToSplitter(storm::storage::sparse::state_type state); + // Initializes the quotient distributions wrt. to the current partition. + void initializeQuotientDistributions(); - // Moves the given state to the position marked by marker1 moves the marker one step further. - void moveStateToMarker1(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock); + // Retrieves whether the given block possibly needs refinement. + bool possiblyNeedsRefinement(bisimulation::Block const& block) const; - // Moves the given state to the position marked by marker2 the marker one step further. - void moveStateToMarker2(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock); + // Splits the given block according to the current quotient distributions. + bool splitBlockAccordingToCurrentQuotientDistributions(bisimulation::Block& block, std::deque*>& splitterQueue); - // Moves the given state to a proper place in the splitter, depending on where the predecessor is located. - void moveStateInSplitter(storm::storage::sparse::state_type predecessor, bisimulation::Block& predecessorBlock, storm::storage::sparse::state_type currentPositionInSplitter, uint_fast64_t& elementsToSkip); + // Retrieves whether the quotient distributions of state 1 are considered to be less than the ones of state 2. + bool quotientDistributionsLess(storm::storage::sparse::state_type state1, storm::storage::sparse::state_type state2); - // Inserts the block into the list of predecessors if it is not already contained. - void insertIntoPredecessorList(bisimulation::Block& predecessorBlock, std::list*>& predecessorBlocks); + // Updates the ordered list of quotient distribution for the given state. + void updateOrderedQuotientDistributions(storm::storage::sparse::state_type state); - // Explores the remaining states of the splitter. - void exploreRemainingStatesOfSplitter(bisimulation::Block& splitter, std::list*>& predecessorBlocks); + // Updates the quotient distributions of the predecessors of the new block by taking the probability mass + // away from the old block. + void updateQuotientDistributionsOfPredecessors(bisimulation::Block const& newBlock, bisimulation::Block const& oldBlock, std::deque*>& splitterQueue); - // Refines the predecessor blocks wrt. strong bisimulation. - void refinePredecessorBlocksOfSplitter(std::list*> const& predecessorBlocks, std::deque*>& splitterQueue); + bool checkQuotientDistributions() const; // A mapping from choice indices to the state state that has this choice. std::vector choiceToStateMapping; - // A vector that holds the probabilities for all nondeterministic choices of all states of going into the - // splitter. This is used by the method that refines a block based on probabilities. - std::vector probabilitiesToCurrentSplitter; + // A vector that holds the quotient distributions for all nondeterministic choices of all states. + std::vector> quotientDistributions; + // A vector that stores for each state the ordered list of quotient distributions. + std::vector const*> orderedQuotientDistributions; }; } } diff --git a/src/utility/ConstantsComparator.cpp b/src/utility/ConstantsComparator.cpp index ed3d82fe4..c3d6018cc 100644 --- a/src/utility/ConstantsComparator.cpp +++ b/src/utility/ConstantsComparator.cpp @@ -35,6 +35,11 @@ namespace storm { return false; } + template + bool ConstantsComparator::isLess(ValueType const& value1, ValueType const& value2) const { + return value1 < value2; + } + ConstantsComparator::ConstantsComparator() : precision(static_cast(storm::settings::generalSettings().getPrecision())) { // Intentionally left empty. } @@ -63,6 +68,10 @@ namespace storm { return value == storm::utility::infinity(); } + bool ConstantsComparator::isLess(float const& value1, float const& value2) const { + return std::abs(value1 - value2) < precision; + } + ConstantsComparator::ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) { // Intentionally left empty. } @@ -90,7 +99,11 @@ namespace storm { bool ConstantsComparator::isConstant(double const& value) const { return true; } - + + bool ConstantsComparator::isLess(double const& value1, double const& value2) const { + return std::abs(value1 - value2) < precision; + } + // Explicit instantiations. template class ConstantsComparator; template class ConstantsComparator; diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h index 7cd03eab3..d30a811d2 100644 --- a/src/utility/ConstantsComparator.h +++ b/src/utility/ConstantsComparator.h @@ -21,6 +21,8 @@ namespace storm { bool isConstant(ValueType const& value) const; bool isInfinity(ValueType const& value) const; + + bool isLess(ValueType const& value1, ValueType const& value2) const; }; // For floats we specialize this class and consider the comparison modulo some predefined precision. @@ -41,6 +43,8 @@ namespace storm { bool isInfinity(float const& value) const; + bool isLess(float const& value1, float const& value2) const; + private: // The precision used for comparisons. float precision; @@ -64,6 +68,8 @@ namespace storm { bool isConstant(double const& value) const; + bool isLess(double const& value1, double const& value2) const; + private: // The precision used for comparisons. double precision; diff --git a/src/utility/storm.h b/src/utility/storm.h index 221545743..d1e32f839 100644 --- a/src/utility/storm.h +++ b/src/utility/storm.h @@ -45,6 +45,7 @@ // Headers for model processing. #include "src/storage/bisimulation/DeterministicModelBisimulationDecomposition.h" +#include "src/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h" // Headers for model checking. #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" @@ -114,7 +115,7 @@ namespace storm { } template - std::shared_ptr performSparseBisimulationMinimization(std::shared_ptr model, std::vector> const& formulas) { + std::shared_ptr performDeterministicSparseBisimulationMinimization(std::shared_ptr model, std::vector> const& formulas) { std::cout << "Performing bisimulation minimization... "; typename storm::storage::DeterministicModelBisimulationDecomposition::Options options; if (!formulas.empty()) { @@ -132,19 +133,40 @@ namespace storm { return model; } + template + std::shared_ptr performNondeterministicSparseBisimulationMinimization(std::shared_ptr model, std::vector> const& formulas) { + std::cout << "Performing bisimulation minimization... "; + typename storm::storage::DeterministicModelBisimulationDecomposition::Options options; + if (!formulas.empty()) { + options = typename storm::storage::NondeterministicModelBisimulationDecomposition::Options(*model, formulas); + } + if (storm::settings::bisimulationSettings().isWeakBisimulationSet()) { + options.type = storm::storage::BisimulationType::Weak; + options.bounded = false; + } + + storm::storage::NondeterministicModelBisimulationDecomposition bisimulationDecomposition(*model, options); + bisimulationDecomposition.computeBisimulationDecomposition(); + model = bisimulationDecomposition.getQuotient(); + std::cout << "done." << std::endl << std::endl; + return model; + } + template, ModelType>::value, bool>::type = 0> std::shared_ptr preprocessModel(std::shared_ptr model, std::vector> const& formulas) { if (storm::settings::generalSettings().isBisimulationSet()) { STORM_LOG_THROW(model->isSparseModel(), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for sparse models."); - STORM_LOG_THROW(model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs and CTMCs."); + STORM_LOG_THROW(model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc) || model->isOfType(storm::models::ModelType::Mdp), storm::exceptions::InvalidSettingsException, "Bisimulation minimization is currently only available for DTMCs, CTMCs and MDPs."); model->reduceToStateBasedRewards(); if (model->isOfType(storm::models::ModelType::Dtmc)) { - return performSparseBisimulationMinimization>(model->template as>(), formulas); + return performDeterministicSparseBisimulationMinimization>(model->template as>(), formulas); + } else if (model->isOfType(storm::models::ModelType::Ctmc)) { + return performDeterministicSparseBisimulationMinimization>(model->template as>(), formulas); } else { - return performSparseBisimulationMinimization>(model->template as>(), formulas); + return performNondeterministicSparseBisimulationMinimization>(model->template as>(), formulas); } }