First version of DFT approximation

Former-commit-id: 1d95ad4914

examples/dft/approx.dft

@@ -0,0 +1,5 @@
+toplevel "A";
+"A" and "B" "C" "D";
+"B" lambda=0.1 dorm=0;
+"C" lambda=10 dorm=0;
+"D" lambda=0.1 dorm=0;

src/builder/ExplicitDFTModelBuilderApprox.cpp

@@ -12,28 +12,29 @@
 namespace storm {
     namespace builder {
 
-        template <typename ValueType, typename StateType>
+        template<typename ValueType, typename StateType>
         ExplicitDFTModelBuilderApprox<ValueType, StateType>::ModelComponents::ModelComponents() : transitionMatrix(), stateLabeling(), markovianStates(), exitRates(), choiceLabeling() {
             // Intentionally left empty.
         }
 
-        template <typename ValueType, typename StateType>
+        template<typename ValueType, typename StateType>
         ExplicitDFTModelBuilderApprox<ValueType, StateType>::ExplicitDFTModelBuilderApprox(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTIndependentSymmetries const& symmetries, bool enableDC) : dft(dft), enableDC(enableDC), stateStorage(((dft.stateVectorSize() / 64) + 1) * 64) {
             // stateVectorSize is bound for size of bitvector
-
             stateGenerationInfo = std::make_shared<storm::storage::DFTStateGenerationInfo>(dft.buildStateGenerationInfo(symmetries));
         }
 
-        template <typename ValueType, typename StateType>
-        std::shared_ptr<storm::models::sparse::Model<ValueType>> ExplicitDFTModelBuilderApprox<ValueType, StateType>::buildModel(LabelOptions const& labelOpts) {
+        template<typename ValueType, typename StateType>
+        void ExplicitDFTModelBuilderApprox<ValueType, StateType>::buildModel(LabelOptions const& labelOpts, double approximationError) {
             STORM_LOG_TRACE("Generating DFT state space");
 
             // Initialize
             StateType currentRowGroup = 0;
             StateType currentRow = 0;
+            size_t pseudoStatesToCheck = 0;
             modelComponents.markovianStates = storm::storage::BitVector(INITIAL_BITVECTOR_SIZE);
             // Create generator
             storm::generator::DftNextStateGenerator<ValueType, StateType> generator(dft, *stateGenerationInfo, enableDC, mergeFailedStates);
+            generator.setApproximationError(approximationError);
             // Create sparse matrix builder
             storm::storage::SparseMatrixBuilder<ValueType> transitionMatrixBuilder(0, 0, 0, false, !generator.isDeterministicModel(), 0);
 
@@ -87,41 +88,55 @@ namespace storm {
                 // Remember that this row group was actually filled with the transitions of a different state
                 setRemapping(currentState->getId(), currentRowGroup);
 
-                // Explore state
-                generator.load(currentState);
-                storm::generator::StateBehavior<ValueType, StateType> behavior = generator.expand(stateToIdCallback);
-
-                STORM_LOG_ASSERT(!behavior.empty(), "Behavior is empty.");
-                setMarkovian(currentRowGroup, behavior.begin()->isMarkovian());
-
                 // If the model is nondeterministic, we need to open a row group.
                 if (!generator.isDeterministicModel()) {
                     transitionMatrixBuilder.newRowGroup(currentRow);
                 }
 
-                // Now add all choices.
-                for (auto const& choice : behavior) {
-                    // Add the probabilistic behavior to the matrix.
-                    for (auto const& stateProbabilityPair : choice) {
-
-                        // Check that pseudo state and its instantiation do not appear together
-                        // TODO Matthias: prove that this is not possible and remove
-                        if (stateProbabilityPair.first >= OFFSET_PSEUDO_STATE) {
-                            StateType newId = stateProbabilityPair.first - OFFSET_PSEUDO_STATE;
-                            STORM_LOG_ASSERT(newId < mPseudoStatesMapping.size(), "Id is not valid.");
-                            if (mPseudoStatesMapping[newId].first > 0) {
-                                // State exists already
-                                newId = mPseudoStatesMapping[newId].first;
-                                for (auto itFind = choice.begin(); itFind != choice.end(); ++itFind) {
-                                    STORM_LOG_ASSERT(itFind->first != newId, "Pseudo state and instantiation occur together in a distribution.");
+                // Try to explore the next state
+                generator.load(currentState);
+
+                if (currentState->isSkip()) {
+                    // Skip the current state
+                    STORM_LOG_TRACE("Skip expansion of state: " << dft.getStateString(currentState));
+                    setMarkovian(currentRowGroup, true);
+                    // Add transition to target state with temporary value 0
+                    // TODO Matthias: what to do when there is no unique target state?
+                    transitionMatrixBuilder.addNextValue(currentRow, failedStateId, storm::utility::zero<ValueType>());
+                    // Remember skipped state
+                    skippedStates[currentRow] = currentState;
+                    ++currentRow;
+                } else {
+                    // Explore the current state
+                    storm::generator::StateBehavior<ValueType, StateType> behavior = generator.expand(stateToIdCallback);
+                    STORM_LOG_ASSERT(!behavior.empty(), "Behavior is empty.");
+                    setMarkovian(currentRowGroup, behavior.begin()->isMarkovian());
+
+                    // Now add all choices.
+                    for (auto const& choice : behavior) {
+                        // Add the probabilistic behavior to the matrix.
+                        for (auto const& stateProbabilityPair : choice) {
+
+                            // Check that pseudo state and its instantiation do not appear together
+                            // TODO Matthias: prove that this is not possible and remove
+                            if (stateProbabilityPair.first >= OFFSET_PSEUDO_STATE) {
+                                StateType newId = stateProbabilityPair.first - OFFSET_PSEUDO_STATE;
+                                STORM_LOG_ASSERT(newId < mPseudoStatesMapping.size(), "Id is not valid.");
+                                if (mPseudoStatesMapping[newId].first > 0) {
+                                    // State exists already
+                                    newId = mPseudoStatesMapping[newId].first;
+                                    for (auto itFind = choice.begin(); itFind != choice.end(); ++itFind) {
+                                        STORM_LOG_ASSERT(itFind->first != newId, "Pseudo state and instantiation occur together in a distribution.");
+                                    }
                                 }
                             }
-                        }
 
-                        transitionMatrixBuilder.addNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second);
+                            transitionMatrixBuilder.addNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second);
+                        }
+                        ++currentRow;
                     }
-                    ++currentRow;
                 }
+
                 ++currentRowGroup;
 
                 if (statesToExplore.empty()) {
@@ -149,6 +164,7 @@ namespace storm {
 
             size_t stateSize = stateStorage.getNumberOfStates() + (mergeFailedStates ? 1 : 0);
             modelComponents.markovianStates.resize(stateSize);
+            modelComponents.deterministicModel = generator.isDeterministicModel();
 
             // Replace pseudo states in matrix
             // TODO Matthias: avoid hack with fixed int type
@@ -169,8 +185,8 @@ namespace storm {
 
             STORM_LOG_TRACE("State remapping: " << stateRemapping);
             STORM_LOG_TRACE("Markovian states: " << modelComponents.markovianStates);
-
             STORM_LOG_DEBUG("Generated " << stateSize << " states");
+            STORM_LOG_DEBUG("Skipped " << skippedStates.size() << " states");
             STORM_LOG_DEBUG("Model is " << (generator.isDeterministicModel() ? "deterministic" : "non-deterministic"));
 
             // Build transition matrix
@@ -222,18 +238,42 @@ namespace storm {
                     }
                 }
             }
+        }
 
-            std::shared_ptr<storm::models::sparse::Model<ValueType>> model;
+        template<typename ValueType, typename StateType>
+        std::shared_ptr<storm::models::sparse::Model<ValueType>> ExplicitDFTModelBuilderApprox<ValueType, StateType>::getModel() {
+            STORM_LOG_ASSERT(skippedStates.size() == 0, "Concrete model has skipped states");
+            return createModel(false);
+        }
+
+        template<typename ValueType, typename StateType>
+        std::shared_ptr<storm::models::sparse::Model<ValueType>> ExplicitDFTModelBuilderApprox<ValueType, StateType>::getModelApproximation(bool lowerBound) {
+            // TODO Matthias: handle case with no skipped states
+            if (lowerBound) {
+                changeMatrixLowerBound(modelComponents.transitionMatrix);
+            } else {
+                changeMatrixUpperBound(modelComponents.transitionMatrix);
+            }
+            return createModel(true);
+        }
 
-            if (generator.isDeterministicModel()) {
+        template<typename ValueType, typename StateType>
+        std::shared_ptr<storm::models::sparse::Model<ValueType>> ExplicitDFTModelBuilderApprox<ValueType, StateType>::createModel(bool copy) {
+            if (modelComponents.deterministicModel) {
                 // Build CTMC
-                model = std::make_shared<storm::models::sparse::Ctmc<ValueType>>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling));
+                if (copy) {
+                    return std::make_shared<storm::models::sparse::Ctmc<ValueType>>(modelComponents.transitionMatrix, modelComponents.stateLabeling);
+
+                } else {
+                    return std::make_shared<storm::models::sparse::Ctmc<ValueType>>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling));
+                }
             } else {
                 // Build MA
                 // Compute exit rates
-                modelComponents.exitRates = std::vector<ValueType>(stateSize);
+                // TODO Matthias: avoid computing multiple times
+                modelComponents.exitRates = std::vector<ValueType>(modelComponents.markovianStates.size());
                 std::vector<typename storm::storage::SparseMatrix<ValueType>::index_type> indices = modelComponents.transitionMatrix.getRowGroupIndices();
-                for (StateType stateIndex = 0; stateIndex < stateSize; ++stateIndex) {
+                for (StateType stateIndex = 0; stateIndex < modelComponents.markovianStates.size(); ++stateIndex) {
                     if (modelComponents.markovianStates[stateIndex]) {
                         modelComponents.exitRates[stateIndex] = modelComponents.transitionMatrix.getRowSum(indices[stateIndex]);
                     } else {
@@ -242,36 +282,52 @@ namespace storm {
                 }
                 STORM_LOG_TRACE("Exit rates: " << modelComponents.exitRates);
 
-                std::shared_ptr<storm::models::sparse::MarkovAutomaton<ValueType>> ma = std::make_shared<storm::models::sparse::MarkovAutomaton<ValueType>>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.markovianStates), std::move(modelComponents.exitRates));
+                std::shared_ptr<storm::models::sparse::MarkovAutomaton<ValueType>> ma;
+                if (copy) {
+                    ma = std::make_shared<storm::models::sparse::MarkovAutomaton<ValueType>>(modelComponents.transitionMatrix, modelComponents.stateLabeling, modelComponents.markovianStates, modelComponents.exitRates);
+                } else {
+                    ma = std::make_shared<storm::models::sparse::MarkovAutomaton<ValueType>>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.markovianStates), std::move(modelComponents.exitRates));
+                }
                 if (ma->hasOnlyTrivialNondeterminism()) {
                     // Markov automaton can be converted into CTMC
-                    model = ma->convertToCTMC();
+                    return ma->convertToCTMC();
                 } else {
-                    model = ma;
+                    return ma;
                 }
             }
-            
-            return model;
-
         }
 
-        template<>
-        bool belowThreshold(double const& number) {
-            return number < 0.1;
+        template<typename ValueType, typename StateType>
+        void ExplicitDFTModelBuilderApprox<ValueType, StateType>::changeMatrixLowerBound(storm::storage::SparseMatrix<ValueType> & matrix) const {
+            // Set lower bound for skipped states
+            for (auto it = skippedStates.begin(); it != skippedStates.end(); ++it) {
+                auto matrixEntry = matrix.getRow(it->first).begin();
+                STORM_LOG_ASSERT(matrixEntry->getColumn() == failedStateId, "Transition has wrong target state.");
+                // Get the lower bound by considering the failure of the BE with the highest rate
+                // The list is sorted by rate, therefore we consider the first BE for the highest failure rate
+                matrixEntry->setValue(it->second->getFailableBERate(0));
+            }
         }
 
-        template<>
-        bool belowThreshold(storm::RationalFunction const& number) {
-            storm::RationalFunction threshold = storm::utility::one<storm::RationalFunction>() / 10;
-            std::cout << number << " < " << threshold << ": " << (number < threshold) << std::endl;
-            std::map<storm::Variable, storm::RationalNumber> mapping;
-            
-            storm::RationalFunction eval(number.evaluate(mapping));
-            std::cout << "Evaluated: " << eval << std::endl;
-            return eval < threshold;
+        template<typename ValueType, typename StateType>
+        void ExplicitDFTModelBuilderApprox<ValueType, StateType>::changeMatrixUpperBound(storm::storage::SparseMatrix<ValueType> & matrix) const {
+            // Set uppper bound for skipped states
+            for (auto it = skippedStates.begin(); it != skippedStates.end(); ++it) {
+                auto matrixEntry = matrix.getRow(it->first).begin();
+                STORM_LOG_ASSERT(matrixEntry->getColumn() == failedStateId, "Transition has wrong target state.");
+                // Get the upper bound by considering the failure of all BEs
+                // The used formula for the rate is 1/( 1/a + 1/b + ...)
+                // TODO Matthias: improve by using closed formula for AND of all BEs
+                ValueType newRate = storm::utility::zero<ValueType>();
+                for (size_t index = 0; index < it->second->nrFailableBEs(); ++index) {
+                    newRate += storm::utility::one<ValueType>() / it->second->getFailableBERate(index);
+                }
+                newRate = storm::utility::one<ValueType>() / newRate;
+                matrixEntry->setValue(newRate);
+            }
         }
 
-        template <typename ValueType, typename StateType>
+        template<typename ValueType, typename StateType>
         StateType ExplicitDFTModelBuilderApprox<ValueType, StateType>::getOrAddStateIndex(DFTStatePointer const& state) {
             StateType stateId;
             bool changed = false;
@@ -329,7 +385,7 @@ namespace storm {
             return stateId;
         }
 
-        template <typename ValueType, typename StateType>
+        template<typename ValueType, typename StateType>
         void ExplicitDFTModelBuilderApprox<ValueType, StateType>::setMarkovian(StateType id, bool markovian) {
             if (id >= modelComponents.markovianStates.size()) {
                 // Resize BitVector
@@ -338,7 +394,7 @@ namespace storm {
             modelComponents.markovianStates.set(id, markovian);
         }
 
-        template <typename ValueType, typename StateType>
+        template<typename ValueType, typename StateType>
         void ExplicitDFTModelBuilderApprox<ValueType, StateType>::setRemapping(StateType id, StateType mappedId) {
             STORM_LOG_ASSERT(id < stateRemapping.size(), "Invalid index for remapping.");
             stateRemapping[id] = mappedId;

src/builder/ExplicitDFTModelBuilderApprox.h

@@ -48,6 +48,9 @@ namespace storm {
 
                 // A vector that stores a labeling for each choice.
                 boost::optional<std::vector<boost::container::flat_set<StateType>>> choiceLabeling;
+
+                // A flag indicating if the model is deterministic.
+                bool deterministicModel;
             };
 
         public:
@@ -68,13 +71,28 @@ namespace storm {
             ExplicitDFTModelBuilderApprox(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTIndependentSymmetries const& symmetries, bool enableDC);
 
             /*!
-             * Build model from Dft.
+             * Build model from DFT.
              *
-             * @param labelOpts Options for labeling.
+             * @param labelOpts          Options for labeling.
+             * @param approximationError Error allowed for approximation.
+             */
+            void buildModel(LabelOptions const& labelOpts, double approximationError = 0.0);
+
+            /*!
+             * Get the built model.
              *
-             * @return Built model (either MA or CTMC).
+             * @return The model built from the DFT.
              */
-            std::shared_ptr<storm::models::sparse::Model<ValueType>> buildModel(LabelOptions const& labelOpts);
+            std::shared_ptr<storm::models::sparse::Model<ValueType>> getModel();
+
+            /*!
+             * Get the built approximation model for either the lower or upper bound.
+             *
+             * @param lowerBound If true, the lower bound model is returned, else the upper bound model
+             *
+             * @return The model built from the DFT.
+             */
+            std::shared_ptr<storm::models::sparse::Model<ValueType>> getModelApproximation(bool lowerBound = true);
 
         private:
 
@@ -103,6 +121,31 @@ namespace storm {
              */
             void setRemapping(StateType id, StateType mappedId);
 
+            /**
+             * Change matrix to reflect the lower approximation bound.
+             *
+             * @param matrix Matrix to change. The change are reflected here.
+             */
+            void changeMatrixLowerBound(storm::storage::SparseMatrix<ValueType> & matrix) const;
+
+            /*!
+             * Change matrix to reflect the upper approximation bound.
+             *
+             * @param matrix Matrix to change. The change are reflected here.
+             */
+            void changeMatrixUpperBound(storm::storage::SparseMatrix<ValueType> & matrix) const;
+
+            /*!
+             * Create the model model from the model components.
+             *
+             * @param copy If true, all elements of the model component are copied (used for approximation). If false
+             *             they are moved to reduce the memory overhead.
+             *
+             * @return The model built from the model components.
+             */
+            std::shared_ptr<storm::models::sparse::Model<ValueType>> createModel(bool copy);
+
+
             // Initial size of the bitvector.
             const size_t INITIAL_BITVECTOR_SIZE = 20000;
             // Offset used for pseudo states.
@@ -118,13 +161,17 @@ namespace storm {
             // Current id for new state
             size_t newIndex = 0;
 
-            //TODO Matthias: remove when everything works
-            std::vector<std::pair<StateType, storm::storage::BitVector>> mPseudoStatesMapping; // vector of (id to concrete state, bitvector)
+            // Mapping from pseudo states to (id of concrete state, bitvector)
+            std::vector<std::pair<StateType, storm::storage::BitVector>> mPseudoStatesMapping;
+
             //TODO Matthias: make changeable
             const bool mergeFailedStates = true;
+
+            // Id of failed state
             size_t failedStateId = 0;
+
+            // Id of initial state
             size_t initialStateIndex = 0;
-            size_t pseudoStatesToCheck = 0;
 
             // Flag indication if dont care propagation should be used.
             bool enableDC = true;
@@ -142,10 +189,10 @@ namespace storm {
             // TODO Matthias: avoid hack with fixed int type
             std::vector<uint_fast64_t> stateRemapping;
 
+            // Holds all skipped states which were not yet expanded. More concrete it is a mapping from matrix indices
+            // to the corresponding skipped state.
+            std::unordered_map<StateType, DFTStatePointer> skippedStates;
         };
-        
-        template<typename ValueType>
-        bool belowThreshold(ValueType const& number);
     }
 }
 

src/generator/DftNextStateGenerator.cpp

@@ -8,7 +8,7 @@ namespace storm {
     namespace generator {
         
         template<typename ValueType, typename StateType>
-        DftNextStateGenerator<ValueType, StateType>::DftNextStateGenerator(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTStateGenerationInfo const& stateGenerationInfo, bool enableDC, bool mergeFailedStates) : mDft(dft), mStateGenerationInfo(stateGenerationInfo), state(nullptr), enableDC(enableDC), mergeFailedStates(mergeFailedStates), comparator() {
+        DftNextStateGenerator<ValueType, StateType>::DftNextStateGenerator(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTStateGenerationInfo const& stateGenerationInfo, bool enableDC, bool mergeFailedStates) : mDft(dft), mStateGenerationInfo(stateGenerationInfo), state(nullptr), enableDC(enableDC), mergeFailedStates(mergeFailedStates) {
             deterministicModel = !mDft.canHaveNondeterminism();
         }
         
@@ -60,7 +60,7 @@ namespace storm {
             Choice<ValueType, StateType> choice(0, !hasDependencies);
 
             // Check for absorbing state
-            if (mDft.hasFailed(state) || mDft.isFailsafe(state) || state->nrFailableBEs() == 0) {
+            if (mDft.hasFailed(state) || mDft.isFailsafe(state) || failableCount == 0) {
                 // Add self loop
                 choice.addProbability(state->getId(), storm::utility::one<ValueType>());
                 STORM_LOG_TRACE("Added self loop for " << state->getId());
@@ -82,30 +82,6 @@ namespace storm {
                 STORM_LOG_ASSERT(nextBEPair.second == hasDependencies, "Failure due to dependencies does not match.");
                 STORM_LOG_TRACE("With the failure of: " << nextBE->name() << " [" << nextBE->id() << "] in " << mDft.getStateString(state));
 
-                /*if (storm::settings::getModule<storm::settings::modules::DFTSettings>().computeApproximation()) {
-                    if (!storm::utility::isZero(exitRate)) {
-                        ValueType rate = nextBE->activeFailureRate();
-                        ValueType div = rate / exitRate;
-                        if (!storm::utility::isZero(exitRate) && belowThreshold(div)) {
-                            // Set transition directly to failed state
-                            auto resultFind = outgoingRates.find(failedIndex);
-                            if (resultFind != outgoingRates.end()) {
-                                // Add to existing transition
-                                resultFind->second += rate;
-                                STORM_LOG_TRACE("Updated transition to " << resultFind->first << " with rate " << rate << " to new rate " << resultFind->second);
-                            } else {
-                                // Insert new transition
-                                outgoingRates.insert(std::make_pair(failedIndex, rate));
-                                STORM_LOG_TRACE("Added transition to " << failedIndex << " with rate " << rate);
-                            }
-                            exitRate += rate;
-                            std::cout << "IGNORE: " << nextBE->name() << " [" << nextBE->id() << "] with rate " << rate << std::endl;
-                            //STORM_LOG_TRACE("Ignore: " << nextBE->name() << " [" << nextBE->id() << "] with rate " << rate);
-                            continue;
-                        }
-                    }
-                }*/
-
                 // Propagate
                 storm::storage::DFTStateSpaceGenerationQueues<ValueType> queues;
 
@@ -139,8 +115,8 @@ namespace storm {
                     continue;
                 }
 
+                // Get the id of the successor state
                 StateType newStateId;
-
                 if (newState->hasFailed(mDft.getTopLevelIndex()) && mergeFailedStates) {
                     // Use unique failed state
                     newStateId = mergeFailedStateId;
@@ -168,17 +144,17 @@ namespace storm {
                 // Set transitions
                 if (hasDependencies) {
                     // Failure is due to dependency -> add non-deterministic choice
-                    std::shared_ptr<storm::storage::DFTDependency<ValueType> const> dependency = mDft.getDependency(state->getDependencyId(currentFailable));
-                    choice.addProbability(newStateId, dependency->probability());
-                    STORM_LOG_TRACE("Added transition to " << newStateId << " with probability " << dependency->probability());
+                    ValueType probability = mDft.getDependency(state->getDependencyId(currentFailable))->probability();
+                    choice.addProbability(newStateId, probability);
+                    STORM_LOG_TRACE("Added transition to " << newStateId << " with probability " << probability);
 
-                    if (!storm::utility::isOne(dependency->probability())) {
+                    if (!storm::utility::isOne(probability)) {
                         // Add transition to state where dependency was unsuccessful
                         DFTStatePointer unsuccessfulState = std::make_shared<storm::storage::DFTState<ValueType>>(*state);
                         unsuccessfulState->letDependencyBeUnsuccessful(currentFailable);
                         // Add state
                         StateType unsuccessfulStateId = stateToIdCallback(unsuccessfulState);
-                        ValueType remainingProbability = storm::utility::one<ValueType>() - dependency->probability();
+                        ValueType remainingProbability = storm::utility::one<ValueType>() - probability;
                         choice.addProbability(unsuccessfulStateId, remainingProbability);
                         STORM_LOG_TRACE("Added transition to " << unsuccessfulStateId << " with remaining probability " << remainingProbability);
                     }
@@ -188,14 +164,21 @@ namespace storm {
                     // Set failure rate according to activation
                     bool isActive = true;
                     if (mDft.hasRepresentant(nextBE->id())) {
-                        // Active must be checked for the state we are coming from as this state is responsible for the
-                        // rate and not the new state we are going to
+                        // Active must be checked for the state we are coming from as this state is responsible for the rate
                         isActive = state->isActive(mDft.getRepresentant(nextBE->id())->id());
                     }
                     ValueType rate = isActive ? nextBE->activeFailureRate() : nextBE->passiveFailureRate();
                     STORM_LOG_ASSERT(!storm::utility::isZero(rate), "Rate is 0.");
                     choice.addProbability(newStateId, rate);
-                    STORM_LOG_TRACE("Added transition to " << newStateId << " with " << (isActive ? "active" : "passive") << " rate " << rate);
+                    STORM_LOG_TRACE("Added transition to " << newStateId << " with " << (isActive ? "active" : "passive") << " failure rate " << rate);
+
+                    // Check if new state needs expansion for approximation
+                    if (approximationError > 0.0) {
+                        if (checkSkipState(newState, rate, choice.getTotalMass(), approximationError)) {
+                            STORM_LOG_TRACE("Will skip state " << newStateId);
+                            newState->markSkip();
+                        }
+                    }
                 }
 
                 ++currentFailable;
@@ -230,6 +213,36 @@ namespace storm {
             return result;
         }
 
+        template<typename ValueType, typename StateType>
+        void DftNextStateGenerator<ValueType, StateType>::setApproximationError(double approximationError) {
+            this->approximationError = approximationError;
+        }
+
+        template<typename ValueType, typename StateType>
+        bool DftNextStateGenerator<ValueType, StateType>::checkSkipState(DFTStatePointer const& state, ValueType rate, ValueType exitRate, double approximationError) {
+            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Approximation works only for double.");
+        }
+
+        template<>
+        bool DftNextStateGenerator<double>::checkSkipState(DFTStatePointer const& state, double rate, double exitRate, double approximationError) {
+            if (mDft.hasFailed(state) || mDft.isFailsafe(state) ||  (state->nrFailableDependencies() == 0 && state->nrFailableBEs() == 0)) {
+                // Do not skip absorbing state
+                return false;
+            }
+            // Skip if the rate to reach this state is low compared to all other outgoing rates from the predecessor
+            return rate/exitRate < approximationError;
+        }
+
+        template<>
+        bool DftNextStateGenerator<storm::RationalFunction>::checkSkipState(DFTStatePointer const& state, storm::RationalFunction rate, storm::RationalFunction exitRate, double approximationError) {
+            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Approximation works only for double.");
+            /*std::cout << (rate / exitRate) << " < " << threshold << ": " << (number < threshold) << std::endl;
+            std::map<storm::Variable, storm::RationalNumber> mapping;
+            storm::RationalFunction eval(number.evaluate(mapping));
+            std::cout << "Evaluated: " << eval << std::endl;
+            return eval < threshold;*/
+        }
+        
         template class DftNextStateGenerator<double>;
         template class DftNextStateGenerator<storm::RationalFunction>;
     }

src/generator/DftNextStateGenerator.h

@@ -42,6 +42,13 @@ namespace storm {
              */
             StateBehavior<ValueType, StateType> createMergeFailedState(StateToIdCallback const& stateToIdCallback);
 
+            /*!
+             * Set a new value for the allowed approximation error.
+             *
+             * @param approximationError Allowed approximation error.
+             */
+            void setApproximationError(double approximationError);
+
         private:
             
             // The dft used for the generation of next states.
@@ -65,8 +72,20 @@ namespace storm {
             // Flag indicating if the model is deterministic.
             bool deterministicModel = true;
 
-            // A comparator used to compare constants.
-            storm::utility::ConstantsComparator<ValueType> comparator;
+            // Allowed approximation error.
+            double approximationError = 0.0;
+
+            /*!
+             * Check if the given state should be skipped for expansion.
+             *
+             * @param state              State to check for expansion
+             * @param rate               Rate of current state
+             * @param exitRate           Exit rates of all outgoing transitions
+             * @param approximationError Allowed approximation error
+             *
+             * @return True, if the given state should be skipped.
+             */
+            bool checkSkipState(DFTStatePointer const& state, ValueType rate, ValueType exitRate, double approximationError);
         };
         
     }

src/modelchecker/dft/DFTModelChecker.cpp

@@ -11,158 +11,213 @@ namespace storm {
 
         template<typename ValueType>
         DFTModelChecker<ValueType>::DFTModelChecker() {
-            // Intentionally left empty.
+            checkResult = storm::utility::zero<ValueType>();
         }
 
         template<typename ValueType>
-        void DFTModelChecker<ValueType>::check(storm::storage::DFT<ValueType> const& origDft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC) {
+        void DFTModelChecker<ValueType>::check(storm::storage::DFT<ValueType> const& origDft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC, double approximationError) {
+            // Initialize
+            this->buildingTime = std::chrono::duration<double>::zero();
+            this->explorationTime = std::chrono::duration<double>::zero();
+            this->bisimulationTime = std::chrono::duration<double>::zero();
+            this->modelCheckingTime = std::chrono::duration<double>::zero();
+            this->totalTime = std::chrono::duration<double>::zero();
+            this->approximationError = approximationError;
             std::chrono::high_resolution_clock::time_point totalStart = std::chrono::high_resolution_clock::now();
+
             // Optimizing DFT
             storm::storage::DFT<ValueType> dft = origDft.optimize();
-            checkResult = checkHelper(dft, formula, symred, allowModularisation, enableDC);
-            totalTime = std::chrono::high_resolution_clock::now() - totalStart;
+
+            // TODO Matthias: check that all paths reach the target state!
+
+            // Checking DFT
+            checkResult = checkHelper(dft, formula, symred, allowModularisation, enableDC, approximationError);
+            this->totalTime = std::chrono::high_resolution_clock::now() - totalStart;
         }
 
         template<typename ValueType>
-        ValueType DFTModelChecker<ValueType>::checkHelper(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC)  {
+        typename DFTModelChecker<ValueType>::dft_result DFTModelChecker<ValueType>::checkHelper(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC, double approximationError)  {
             STORM_LOG_TRACE("Check helper called");
-            bool invResults = false;
-            std::vector<storm::storage::DFT<ValueType>> dfts = {dft};
-            std::vector<ValueType> res;
-            size_t nrK = 0; // K out of M
-            size_t nrM = 0; // K out of M
+            bool modularisationPossible = allowModularisation;
 
             // Try modularisation
-            if(allowModularisation) {
-                if(dft.topLevelType() == storm::storage::DFTElementType::AND) {
-                    STORM_LOG_TRACE("top modularisation called AND");
-                    dfts = dft.topModularisation();
-                    STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
-                    nrK = dfts.size();
-                    nrM = dfts.size();
-                }
-                if(dft.topLevelType() == storm::storage::DFTElementType::OR) {
-                    STORM_LOG_TRACE("top modularisation called OR");
-                    dfts = dft.topModularisation();
-                    STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
-                    nrK = 0;
-                    nrM = dfts.size();
-                    invResults = true;
-                }
-                if(dft.topLevelType() == storm::storage::DFTElementType::VOT) {
-                    STORM_LOG_TRACE("top modularisation called VOT");
-                    dfts = dft.topModularisation();
-                    STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
-                    nrK = std::static_pointer_cast<storm::storage::DFTVot<ValueType> const>(dft.getTopLevelGate())->threshold();
-                    nrM = dfts.size();
-                    if(nrK <= nrM/2) {
-                        nrK -= 1;
+            if(modularisationPossible) {
+                bool invResults = false;
+                std::vector<storm::storage::DFT<ValueType>> dfts;
+                size_t nrK = 0; // K out of M
+                size_t nrM = 0; // K out of M
+
+                switch (dft.topLevelType()) {
+                    case storm::storage::DFTElementType::AND:
+                        STORM_LOG_TRACE("top modularisation called AND");
+                        dfts = dft.topModularisation();
+                        STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
+                        nrK = dfts.size();
+                        nrM = dfts.size();
+                        modularisationPossible = dfts.size() > 1;
+                        break;
+                    case storm::storage::DFTElementType::OR:
+                        STORM_LOG_TRACE("top modularisation called OR");
+                        dfts = dft.topModularisation();
+                        STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
+                        nrK = 0;
+                        nrM = dfts.size();
                         invResults = true;
-                    }
+                        modularisationPossible = dfts.size() > 1;
+                        break;
+                    case storm::storage::DFTElementType::VOT:
+                        STORM_LOG_TRACE("top modularisation called VOT");
+                        dfts = dft.topModularisation();
+                        STORM_LOG_TRACE("Modularsation into " << dfts.size() << " submodules.");
+                        nrK = std::static_pointer_cast<storm::storage::DFTVot<ValueType> const>(dft.getTopLevelGate())->threshold();
+                        nrM = dfts.size();
+                        if(nrK <= nrM/2) {
+                            nrK -= 1;
+                            invResults = true;
+                        }
+                        modularisationPossible = dfts.size() > 1;
+                        break;
+                    default:
+                        // No static gate -> no modularisation applicable
+                        modularisationPossible = false;
+                        break;
                 }
-                if(dfts.size() > 1) {
+
+                if(modularisationPossible) {
                     STORM_LOG_TRACE("Recursive CHECK Call");
+                    // TODO Matthias: enable modularisation for approximation
+                    STORM_LOG_ASSERT(approximationError == 0.0, "Modularisation not possible for approximation.");
+
+                    // Recursively call model checking
+                    std::vector<ValueType> res;
                     for(auto const ft : dfts) {
-                        res.push_back(checkHelper(ft, formula, symred, allowModularisation, enableDC));
+                        dft_result ftResult = checkHelper(ft, formula, symred, true, enableDC, 0.0);
+                        res.push_back(boost::get<ValueType>(ftResult));
                     }
-                }
-            }
 
-            if(res.empty()) {
-                // Model based modularisation.
-                auto const& models = buildMarkovModels(dfts, formula, symred, enableDC);
-                
-                for (auto const& model : models) {
-                    // Model checking
-                    STORM_LOG_INFO("Model checking...");
-                    std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now();
-                    std::unique_ptr<storm::modelchecker::CheckResult> result(storm::verifySparseModel(model, formula));
-                    STORM_LOG_INFO("Model checking done.");
-                    STORM_LOG_ASSERT(result, "Result does not exist.");
-                    result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
-                    modelCheckingTime += std::chrono::high_resolution_clock::now() -  modelCheckingStart;
-                    res.push_back(result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second);
+                    // Combine modularisation results
+                    STORM_LOG_TRACE("Combining all results... K=" << nrK << "; M=" << nrM << "; invResults=" << (invResults?"On":"Off"));
+                    ValueType result = storm::utility::zero<ValueType>();
+                    int limK = invResults ? -1 : nrM+1;
+                    int chK = invResults ? -1 : 1;
+                    for(int cK = nrK; cK != limK; cK += chK ) {
+                        STORM_LOG_ASSERT(cK >= 0, "ck negative.");
+                        size_t permutation = smallestIntWithNBitsSet(static_cast<size_t>(cK));
+                        do {
+                            STORM_LOG_TRACE("Permutation="<<permutation);
+                            ValueType permResult = storm::utility::one<ValueType>();
+                            for(size_t i = 0; i < res.size(); ++i) {
+                                if(permutation & (1 << i)) {
+                                    permResult *= res[i];
+                                } else {
+                                    permResult *= storm::utility::one<ValueType>() - res[i];
+                                }
+                            }
+                            STORM_LOG_TRACE("Result for permutation:"<<permResult);
+                            permutation = nextBitPermutation(permutation);
+                            result += permResult;
+                        } while(permutation < (1 << nrM) && permutation != 0);
+                    }
+                    if(invResults) {
+                        result = storm::utility::one<ValueType>() - result;
+                    }
+                    return result;
                 }
             }
-            if(nrM <= 1) {
-                // No modularisation done.
-                STORM_LOG_ASSERT(res.size()==1, "Result size not 1.");
-                return res[0];
-            }
 
-            STORM_LOG_TRACE("Combining all results... K=" << nrK << "; M=" << nrM << "; invResults=" << (invResults?"On":"Off"));
-            ValueType result = storm::utility::zero<ValueType>();
-            int limK = invResults ? -1 : nrM+1;
-            int chK = invResults ? -1 : 1;
-            for(int cK = nrK; cK != limK; cK += chK ) {
-                STORM_LOG_ASSERT(cK >= 0, "ck negative.");
-                size_t permutation = smallestIntWithNBitsSet(static_cast<size_t>(cK));
-                do {
-                    STORM_LOG_TRACE("Permutation="<<permutation);
-                    ValueType permResult = storm::utility::one<ValueType>();
-                    for(size_t i = 0; i < res.size(); ++i) {
-                        if(permutation & (1 << i)) {
-                            permResult *= res[i];
-                        } else {
-                            permResult *= storm::utility::one<ValueType>() - res[i];
-                        }
-                    }
-                    STORM_LOG_TRACE("Result for permutation:"<<permResult);
-                    permutation = nextBitPermutation(permutation);
-                    result += permResult;
-                } while(permutation < (1 << nrM) && permutation != 0);
-            }
-            if(invResults) {
-                return storm::utility::one<ValueType>() - result;
-            }
-            return result;
+            // If we are here, no modularisation was possible
+            STORM_LOG_ASSERT(!modularisationPossible, "Modularisation should not be possible.");
+            return checkDFT(dft, formula, symred, enableDC, approximationError);
         }
 
         template<typename ValueType>
-        std::vector<ValueType> DFTModelChecker<ValueType>::checkModel(std::vector<storm::storage::DFT<ValueType>> const& dfts, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError) {
-            std::vector<ValueType> results;
+        typename DFTModelChecker<ValueType>::dft_result DFTModelChecker<ValueType>::checkDFT(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError) {
+            std::chrono::high_resolution_clock::time_point buildingStart = std::chrono::high_resolution_clock::now();
 
-            auto const& models = buildMarkovModels(dfts, formula, symred, enableDC);
-
-            for (auto const& model : models) {
-                // Model checking
-                STORM_LOG_INFO("Model checking...");
-                std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now();
-                std::unique_ptr<storm::modelchecker::CheckResult> result(storm::verifySparseModel(model, formula));
-                STORM_LOG_INFO("Model checking done.");
-                STORM_LOG_ASSERT(result, "Result does not exist.");
-                result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
-                modelCheckingTime += std::chrono::high_resolution_clock::now() -  modelCheckingStart;
-                results.push_back(result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second);
+            // Find symmetries
+            std::map<size_t, std::vector<std::vector<size_t>>> emptySymmetry;
+            storm::storage::DFTIndependentSymmetries symmetries(emptySymmetry);
+            if(symred) {
+                auto colouring = dft.colourDFT();
+                symmetries = dft.findSymmetries(colouring);
+                STORM_LOG_INFO("Found " << symmetries.groups.size() << " symmetries.");
+                STORM_LOG_TRACE("Symmetries: " << std::endl << symmetries);
             }
-            return results;
-        }
+            std::chrono::high_resolution_clock::time_point buildingEnd = std::chrono::high_resolution_clock::now();
+            buildingTime += buildingEnd - buildingStart;
 
-        template<typename ValueType>
-        std::vector<std::shared_ptr<storm::models::sparse::Model<ValueType>>> DFTModelChecker<ValueType>::buildMarkovModels(std::vector<storm::storage::DFT<ValueType>> const& dfts, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError) {
-            std::vector<std::shared_ptr<storm::models::sparse::Model<ValueType>>> models;
-            for(auto& dft : dfts) {
-                std::chrono::high_resolution_clock::time_point buildingStart = std::chrono::high_resolution_clock::now();
-
-                std::map<size_t, std::vector<std::vector<size_t>>> emptySymmetry;
-                storm::storage::DFTIndependentSymmetries symmetries(emptySymmetry);
-                if(symred) {
-                    auto colouring = dft.colourDFT();
-                    symmetries = dft.findSymmetries(colouring);
-                    STORM_LOG_INFO("Found " << symmetries.groups.size() << " symmetries.");
-                    STORM_LOG_TRACE("Symmetries: " << std::endl << symmetries);
-                }
-                std::chrono::high_resolution_clock::time_point buildingEnd = std::chrono::high_resolution_clock::now();
-                buildingTime += buildingEnd - buildingStart;
+            if (approximationError > 0.0) {
+                // Build approximate Markov Automata for lower and upper bound
+                double currentApproximationError = approximationError;
+                approximation_result approxResult = std::make_pair(storm::utility::zero<ValueType>(), storm::utility::zero<ValueType>());
+                std::chrono::high_resolution_clock::time_point explorationStart;
+                std::shared_ptr<storm::models::sparse::Model<ValueType>> model;
+                storm::builder::ExplicitDFTModelBuilderApprox<ValueType> builder(dft, symmetries, enableDC);
+                typename storm::builder::ExplicitDFTModelBuilderApprox<ValueType>::LabelOptions labeloptions; // TODO initialize this with the formula
 
-                // Building Markov Automaton
+                size_t iteration = 0;
+                do {
+                    // Iteratively build finer models
+                    // TODO Matthias: implement refinement
+                    STORM_LOG_ASSERT(iteration < 1, "Iterative refinement not yet implemented.");
+                    explorationStart = std::chrono::high_resolution_clock::now();
+                    STORM_LOG_INFO("Building model...");
+                    // TODO Matthias refine model using existing model and MC results
+                    currentApproximationError = pow(0.1, iteration) * approximationError;
+                    builder.buildModel(labeloptions, currentApproximationError);
+                    // TODO Matthias: possible to do bisimulation on approximated model and not on concrete one?
+
+                    // Build model for lower bound
+                    STORM_LOG_INFO("Getting model for lower bound...");
+                    model = builder.getModelApproximation(true);
+                    //model->printModelInformationToStream(std::cout);
+                    STORM_LOG_INFO("No. states (Explored): " << model->getNumberOfStates());
+                    STORM_LOG_INFO("No. transitions (Explored): " << model->getNumberOfTransitions());
+                    if (model->getNumberOfStates() <= 15) {
+                        STORM_LOG_TRACE("Transition matrix: " << std::endl << model->getTransitionMatrix());
+                    } else {
+                        STORM_LOG_TRACE("Transition matrix: too big to print");
+                    }
+                    explorationTime += std::chrono::high_resolution_clock::now() - explorationStart;
+                    // Check lower bound
+                    std::unique_ptr<storm::modelchecker::CheckResult> result = checkModel(model, formula);
+                    result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
+                    approxResult.first = result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second;
+
+                    // Build model for upper bound
+                    STORM_LOG_INFO("Getting model for upper bound...");
+                    explorationStart = std::chrono::high_resolution_clock::now();
+                    model = builder.getModelApproximation(false);
+                    //model->printModelInformationToStream(std::cout);
+                    STORM_LOG_INFO("No. states (Explored): " << model->getNumberOfStates());
+                    STORM_LOG_INFO("No. transitions (Explored): " << model->getNumberOfTransitions());
+                    if (model->getNumberOfStates() <= 15) {
+                        STORM_LOG_TRACE("Transition matrix: " << std::endl << model->getTransitionMatrix());
+                    } else {
+                        STORM_LOG_TRACE("Transition matrix: too big to print");
+                    }
+                    explorationTime += std::chrono::high_resolution_clock::now() - explorationStart;
+                    // Check upper bound
+                    result = checkModel(model, formula);
+                    result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
+                    approxResult.second = result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second;
+
+                    ++iteration;
+                    STORM_LOG_TRACE("Result after iteration " << iteration << ": (" << approxResult.first << ", " << approxResult.second << ")");
+                } while (!isApproximationSufficient(approxResult.first, approxResult.second, approximationError));
+
+                STORM_LOG_INFO("Finished approximation after " << iteration << " iteration" << (iteration > 1 ? "s." : "."));
+                return approxResult;
+            } else {
+                // Build a single Markov Automaton
                 STORM_LOG_INFO("Building Model...");
                 std::shared_ptr<storm::models::sparse::Model<ValueType>> model;
                 // TODO Matthias: use only one builder if everything works again
-                if (storm::settings::getModule<storm::settings::modules::DFTSettings>().computeApproximation()) {
+                if (approximationError >= 0.0) {
                     storm::builder::ExplicitDFTModelBuilderApprox<ValueType> builder(dft, symmetries, enableDC);
                     typename storm::builder::ExplicitDFTModelBuilderApprox<ValueType>::LabelOptions labeloptions; // TODO initialize this with the formula
-                    model = builder.buildModel(labeloptions);
+                    builder.buildModel(labeloptions);
+                    model = builder.getModel();
                 } else {
                     storm::builder::ExplicitDFTModelBuilder<ValueType> builder(dft, symmetries, enableDC);
                     typename storm::builder::ExplicitDFTModelBuilder<ValueType>::LabelOptions labeloptions; // TODO initialize this with the formula
@@ -171,21 +226,45 @@ namespace storm {
                 //model->printModelInformationToStream(std::cout);
                 STORM_LOG_INFO("No. states (Explored): " << model->getNumberOfStates());
                 STORM_LOG_INFO("No. transitions (Explored): " << model->getNumberOfTransitions());
-                std::chrono::high_resolution_clock::time_point explorationEnd = std::chrono::high_resolution_clock::now();
-                explorationTime += explorationEnd -buildingEnd;
+                explorationTime += std::chrono::high_resolution_clock::now() - buildingEnd;
 
-                // Bisimulation
-                if (model->isOfType(storm::models::ModelType::Ctmc) && storm::settings::getModule<storm::settings::modules::GeneralSettings>().isBisimulationSet()) {
-                    STORM_LOG_INFO("Bisimulation...");
-                    model =  storm::performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), {formula}, storm::storage::BisimulationType::Weak)->template as<storm::models::sparse::Ctmc<ValueType>>();
-                    //model->printModelInformationToStream(std::cout);
-                }
+                // Model checking
+                std::unique_ptr<storm::modelchecker::CheckResult> result = checkModel(model, formula);
+                result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
+                return result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second;
+            }
+        }
+
+        template<typename ValueType>
+        std::unique_ptr<storm::modelchecker::CheckResult> DFTModelChecker<ValueType>::checkModel(std::shared_ptr<storm::models::sparse::Model<ValueType>>& model, std::shared_ptr<const storm::logic::Formula> const& formula) {
+            // Bisimulation
+            std::chrono::high_resolution_clock::time_point bisimulationStart = std::chrono::high_resolution_clock::now();
+            if (model->isOfType(storm::models::ModelType::Ctmc) && storm::settings::getModule<storm::settings::modules::GeneralSettings>().isBisimulationSet()) {
+                STORM_LOG_INFO("Bisimulation...");
+                model =  storm::performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), {formula}, storm::storage::BisimulationType::Weak)->template as<storm::models::sparse::Ctmc<ValueType>>();
                 STORM_LOG_INFO("No. states (Bisimulation): " << model->getNumberOfStates());
                 STORM_LOG_INFO("No. transitions (Bisimulation): " << model->getNumberOfTransitions());
-                bisimulationTime += std::chrono::high_resolution_clock::now() - explorationEnd;
-                models.push_back(model);
             }
-            return models;
+            std::chrono::high_resolution_clock::time_point bisimulationEnd = std::chrono::high_resolution_clock::now();
+            bisimulationTime += bisimulationEnd - bisimulationStart;
+
+            // Check the model
+            STORM_LOG_INFO("Model checking...");
+            std::unique_ptr<storm::modelchecker::CheckResult> result(storm::verifySparseModel(model, formula));
+            STORM_LOG_INFO("Model checking done.");
+            STORM_LOG_ASSERT(result, "Result does not exist.");
+            modelCheckingTime += std::chrono::high_resolution_clock::now() - bisimulationEnd;
+            return result;
+        }
+
+        template<typename ValueType>
+        bool DFTModelChecker<ValueType>::isApproximationSufficient(ValueType lowerBound, ValueType upperBound, double approximationError) {
+            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Approximation works only for double.");
+        }
+
+        template<>
+        bool DFTModelChecker<double>::isApproximationSufficient(double lowerBound, double upperBound, double approximationError) {
+            return upperBound - lowerBound <= approximationError;
         }
 
         template<typename ValueType>
@@ -201,7 +280,13 @@ namespace storm {
         template<typename ValueType>
         void DFTModelChecker<ValueType>::printResult(std::ostream& os) {
             os << "Result: [";
-            os << checkResult << "]" << std::endl;
+            if (this->approximationError > 0.0) {
+                approximation_result result = boost::get<approximation_result>(checkResult);
+                os << "(" << result.first << ", " << result.second << ")";
+            } else {
+                os << boost::get<ValueType>(checkResult);
+            }
+            os << "]" << std::endl;
         }
 
 

src/modelchecker/dft/DFTModelChecker.h

@@ -11,20 +11,23 @@
 namespace storm {
     namespace modelchecker {
 
-        /**
+        /*!
          * Analyser for DFTs.
          */
         template<typename ValueType>
         class DFTModelChecker {
 
+            typedef std::pair<ValueType, ValueType> approximation_result;
+            typedef boost::variant<ValueType, approximation_result> dft_result;
+
         public:
 
-            /**
+            /*!
              * Constructor.
              */
             DFTModelChecker();
 
-            /**
+            /*!
              * Main method for checking DFTs.
              *
              * @param origDft             Original DFT
@@ -32,17 +35,18 @@ namespace storm {
              * @param symred              Flag indicating if symmetry reduction should be used
              * @param allowModularisation Flag indication if modularisation is allowed
              * @param enableDC            Flag indicating if dont care propagation should be used
+             * @param approximationError  Error allowed for approximation. Value 0 indicates no approximation
              */
-            void check(storm::storage::DFT<ValueType> const& origDft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred = true, bool allowModularisation = true, bool enableDC = true);
+            void check(storm::storage::DFT<ValueType> const& origDft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred = true, bool allowModularisation = true, bool enableDC = true, double approximationError = 0.0);
 
-            /**
+            /*!
              * Print timings of all operations to stream.
              *
              * @param os Output stream to write to.
              */
             void printTimings(std::ostream& os = std::cout);
 
-            /**
+            /*!
              * Print result to stream.
              *
              * @param os Output stream to write to.
@@ -50,16 +54,21 @@ namespace storm {
             void printResult(std::ostream& os = std::cout);
 
         private:
+
             // Timing values
             std::chrono::duration<double> buildingTime = std::chrono::duration<double>::zero();
             std::chrono::duration<double> explorationTime = std::chrono::duration<double>::zero();
             std::chrono::duration<double> bisimulationTime = std::chrono::duration<double>::zero();
             std::chrono::duration<double> modelCheckingTime = std::chrono::duration<double>::zero();
             std::chrono::duration<double> totalTime = std::chrono::duration<double>::zero();
+
             // Model checking result
-            ValueType checkResult = storm::utility::zero<ValueType>();
+            dft_result checkResult;
+
+            // Allowed error bound for approximation
+            double approximationError;
 
-            /**
+            /*!
              * Internal helper for model checking a DFT.
              *
              * @param dft                 DFT
@@ -67,36 +76,45 @@ namespace storm {
              * @param symred              Flag indicating if symmetry reduction should be used
              * @param allowModularisation Flag indication if modularisation is allowed
              * @param enableDC            Flag indicating if dont care propagation should be used
+             * @param approximationError  Error allowed for approximation. Value 0 indicates no approximation
              *
-             * @return Model checking result
+             * @return Model checking result (or in case of approximation two results for lower and upper bound)
              */
-            ValueType checkHelper(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC);
+            dft_result checkHelper(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool allowModularisation, bool enableDC, double approximationError);
 
-            /**
-             * Check the Dfts via model checking.
+            /*!
+             * Check model generated from DFT.
              *
-             * @param dfts               Vector of Dfts
+             * @param dft                The DFT
              * @param formula            Formula to check for
              * @param symred             Flag indicating if symmetry reduction should be used
              * @param enableDC           Flag indicating if dont care propagation should be used
              * @param approximationError Error allowed for approximation. Value 0 indicates no approximation
              *
-             * @return Vector of results for each model
+             * @return Model checking result
              */
-            std::vector<ValueType> checkModel(std::vector<storm::storage::DFT<ValueType>> const& dfts, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError = 0.0);
+            dft_result checkDFT(storm::storage::DFT<ValueType> const& dft, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError = 0.0);
 
-            /**
-             * Build markov models from DFTs.
+            /*!
+             * Check the given markov model for the given property.
              *
-             * @param dfts               Vector of Dfts
-             * @param formula            Formula to check for
-             * @param symred             Flag indicating if symmetry reduction should be used
-             * @param enableDC           Flag indicating if dont care propagation should be used
-             * @param approximationError Error allowed for approximation. Value 0 indicates no approximation
+             * @param model   Model to check
+             * @param formula Formula to check for
+             *
+             * @return Model checking result
+             */
+            std::unique_ptr<storm::modelchecker::CheckResult> checkModel(std::shared_ptr<storm::models::sparse::Model<ValueType>>& model, std::shared_ptr<const storm::logic::Formula> const& formula);
+
+            /*!
+             * Checks if the computed approximation is sufficient, i.e. upperBound - lowerBound <= approximationError.
+             *
+             * @param lowerBound         The lower bound on the result.
+             * @param upperBound         The upper bound on the result.
+             * @param approximationError The allowed error for approximating.
              *
-             * @return Vector of markov models corresponding to DFTs.
+             * @return True, if the approximation is sufficient.
              */
-            std::vector<std::shared_ptr<storm::models::sparse::Model<ValueType>>> buildMarkovModels(std::vector<storm::storage::DFT<ValueType>> const& dfts, std::shared_ptr<const storm::logic::Formula> const& formula, bool symred, bool enableDC, double approximationError = 0.0);
+            bool isApproximationSufficient(ValueType lowerBound, ValueType upperBound, double approximationError);
 
         };
     }

src/settings/modules/DFTSettings.cpp

@@ -21,8 +21,8 @@ namespace storm {
             const std::string DFTSettings::symmetryReductionOptionShortName = "symred";
             const std::string DFTSettings::modularisationOptionName = "modularisation";
             const std::string DFTSettings::disableDCOptionName = "disabledc";
-            const std::string DFTSettings::computeApproximationOptionName = "approximation";
-            const std::string DFTSettings::computeApproximationOptionShortName = "approx";
+            const std::string DFTSettings::approximationErrorOptionName = "approximation";
+            const std::string DFTSettings::approximationErrorOptionShortName = "approx";
             const std::string DFTSettings::propExpectedTimeOptionName = "expectedtime";
             const std::string DFTSettings::propExpectedTimeOptionShortName = "mttf";
             const std::string DFTSettings::propProbabilityOptionName = "probability";
@@ -39,7 +39,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, symmetryReductionOptionName, false, "Exploit symmetric structure of model.").setShortName(symmetryReductionOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, modularisationOptionName, false, "Use modularisation (not applicable for expected time).").build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, disableDCOptionName, false, "Disable Dont Care propagation.").build());
-                this->addOption(storm::settings::OptionBuilder(moduleName, computeApproximationOptionName, false, "Compute an approximation.").setShortName(computeApproximationOptionShortName).build());
+                this->addOption(storm::settings::OptionBuilder(moduleName, approximationErrorOptionName, false, "Approximation error allowed.").setShortName(approximationErrorOptionShortName).addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("error", "The approximation error to use.").addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleGreaterValidatorIncluding(0.0)).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, propExpectedTimeOptionName, false, "Compute expected time of system failure.").setShortName(propExpectedTimeOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, propProbabilityOptionName, false, "Compute probability of system failure.").build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, propTimeBoundOptionName, false, "Compute probability of system failure up to given timebound.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("time", "The timebound to use.").addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleGreaterValidatorExcluding(0.0)).build()).build());
@@ -70,8 +70,12 @@ namespace storm {
                 return this->getOption(disableDCOptionName).getHasOptionBeenSet();
             }
 
-            bool DFTSettings::computeApproximation() const {
-                return this->getOption(computeApproximationOptionName).getHasOptionBeenSet();
+            bool DFTSettings::isApproximationErrorSet() const {
+                return this->getOption(approximationErrorOptionName).getHasOptionBeenSet();
+            }
+
+            double DFTSettings::getApproximationError() const {
+                return this->getOption(approximationErrorOptionName).getArgumentByName("error").getValueAsDouble();
             }
 
             bool DFTSettings::usePropExpectedTime() const {

src/settings/modules/DFTSettings.h

@@ -33,9 +33,6 @@ namespace storm {
                  */
                 std::string getDftFilename() const;
 
-                //expectedtime, probability, timebound, prob
-                //min, max
-                
                 /*!
                  * Retrieves whether the option to use symmetry reduction is set.
                  *
@@ -62,7 +59,14 @@ namespace storm {
                  *
                  * @return True iff the option was set.
                  */
-                bool computeApproximation() const;
+                bool isApproximationErrorSet() const;
+
+                /*!
+                 * Retrieves the error allowed for approximation the model checking result.
+                 *
+                 * @return The allowed errorbound.
+                 */
+                double getApproximationError() const;
 
                 /*!
                  * Retrieves whether the property expected time should be used.
@@ -129,8 +133,8 @@ namespace storm {
                 static const std::string symmetryReductionOptionShortName;
                 static const std::string modularisationOptionName;
                 static const std::string disableDCOptionName;
-                static const std::string computeApproximationOptionName;
-                static const std::string computeApproximationOptionShortName;
+                static const std::string approximationErrorOptionName;
+                static const std::string approximationErrorOptionShortName;
                 static const std::string propExpectedTimeOptionName;
                 static const std::string propExpectedTimeOptionShortName;
                 static const std::string propProbabilityOptionName;

src/storage/dft/DFTState.cpp

@@ -204,8 +204,14 @@ namespace storm {
         }
 
         template<typename ValueType>
-        std::pair<std::shared_ptr<DFTBE<ValueType> const>, bool> DFTState<ValueType>::letNextBEFail(size_t index)
-        {
+        ValueType DFTState<ValueType>::getFailableBERate(size_t index) const {
+            STORM_LOG_ASSERT(index < nrFailableBEs(), "Index invalid.");
+            // TODO Matthias: get passive failure rate when applicable
+            return mDft.getBasicElement(mIsCurrentlyFailableBE[index])->activeFailureRate();
+        }
+
+        template<typename ValueType>
+        std::pair<std::shared_ptr<DFTBE<ValueType> const>, bool> DFTState<ValueType>::letNextBEFail(size_t index) {
             STORM_LOG_TRACE("currently failable: " << getCurrentlyFailableString());
             if (nrFailableDependencies() > 0) {
                 // Consider failure due to dependency

src/storage/dft/DFTState.h

@@ -31,6 +31,7 @@ namespace storm {
             bool mValid = true;
             const DFT<ValueType>& mDft;
             const DFTStateGenerationInfo& mStateGenerationInfo;
+            bool mSkip = false;
             
         public:
             DFTState(DFT<ValueType> const& dft, DFTStateGenerationInfo const& stateGenerationInfo, size_t id);
@@ -96,6 +97,14 @@ namespace storm {
                 return !mValid;
             }
 
+            void markSkip() {
+                mSkip = true;
+            }
+
+            bool isSkip() const {
+                return mSkip;
+            }
+            
             storm::storage::BitVector const& status() const {
                 return mStatus;
             }
@@ -145,6 +154,14 @@ namespace storm {
                 return mIsCurrentlyFailableBE.size();
             }
 
+            /** Get the failure rate of the currently failable BE on the given index.
+             *
+             * @param index Index of BE in list of currently failable BEs.
+             *
+             * @return Failure rate of the BE.
+             */
+            ValueType getFailableBERate(size_t index) const;
+
             size_t nrFailableDependencies() const {
                 return mFailableDependencies.size();
             }

src/storm-dyftee.cpp

@@ -33,7 +33,7 @@
  * @param property PCTC formula capturing the property to check.
  */
 template <typename ValueType>
-void analyzeDFT(std::string filename, std::string property, bool symred = false, bool allowModularisation = false, bool enableDC = true) {
+void analyzeDFT(std::string filename, std::string property, bool symred, bool allowModularisation, bool enableDC, double approximationError) {
     std::cout << "Running DFT analysis on file " << filename << " with property " << property << std::endl;
 
     storm::parser::DFTGalileoParser<ValueType> parser;
@@ -42,11 +42,16 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false,
     STORM_LOG_ASSERT(formulas.size() == 1, "Wrong number of formulas.");
 
     storm::modelchecker::DFTModelChecker<ValueType> modelChecker;
-    modelChecker.check(dft, formulas[0], symred, allowModularisation, enableDC);
+    modelChecker.check(dft, formulas[0], symred, allowModularisation, enableDC, approximationError);
     modelChecker.printTimings();
     modelChecker.printResult();
 }
 
+/*!
+ * Analyze the DFT with use of SMT solving.
+ *
+ * @param filename Path to DFT file in Galileo format.
+ */
 template<typename ValueType>
 void analyzeWithSMT(std::string filename) {
     std::cout << "Running DFT analysis on file " << filename << " with use of SMT" << std::endl;
@@ -163,11 +168,16 @@ int main(const int argc, const char** argv) {
         
         STORM_LOG_ASSERT(!pctlFormula.empty(), "Pctl formula empty.");
 
+        double approximationError = 0.0;
+        if (dftSettings.isApproximationErrorSet()) {
+            approximationError = dftSettings.getApproximationError();
+        }
+
         // From this point on we are ready to carry out the actual computations.
         if (parametric) {
-            analyzeDFT<storm::RationalFunction>(dftSettings.getDftFilename(), pctlFormula, dftSettings.useSymmetryReduction(), allowModular && dftSettings.useModularisation(), !dftSettings.isDisableDC() );
+            analyzeDFT<storm::RationalFunction>(dftSettings.getDftFilename(), pctlFormula, dftSettings.useSymmetryReduction(), allowModular && dftSettings.useModularisation(), !dftSettings.isDisableDC(), approximationError);
         } else {
-            analyzeDFT<double>(dftSettings.getDftFilename(), pctlFormula, dftSettings.useSymmetryReduction(), allowModular && dftSettings.useModularisation(), !dftSettings.isDisableDC());
+            analyzeDFT<double>(dftSettings.getDftFilename(), pctlFormula, dftSettings.useSymmetryReduction(), allowModular && dftSettings.useModularisation(), !dftSettings.isDisableDC(), approximationError);
         }
         
         // All operations have now been performed, so we clean up everything and terminate.