Better upper bound for independent subtrees

Former-commit-id: 64f5a1ca60

src/builder/ExplicitDFTModelBuilderApprox.cpp

@@ -39,6 +39,61 @@ namespace storm {
             // Intentionally left empty.
             // TODO Matthias: remove again
             heuristic = storm::builder::ApproximationHeuristic::PROBABILITY;
+
+            // Compute independent subtrees
+            if (dft.topLevelType() == storm::storage::DFTElementType::OR) {
+                // We only support this for approximation with top level element OR
+                for (auto const& child : dft.getGate(dft.getTopLevelIndex())->children()) {
+                    // Consider all children of the top level gate
+                    std::vector<size_t> isubdft;
+                    if (child->nrParents() > 1 || child->hasOutgoingDependencies()) {
+                        STORM_LOG_TRACE("child " << child->name() << "does not allow modularisation.");
+                        isubdft.clear();
+                    } else if (dft.isGate(child->id())) {
+                        isubdft = dft.getGate(child->id())->independentSubDft(false);
+                    } else {
+                        STORM_LOG_ASSERT(dft.isBasicElement(child->id()), "Child is no BE.");
+                        if(dft.getBasicElement(child->id())->hasIngoingDependencies()) {
+                            STORM_LOG_TRACE("child " << child->name() << "does not allow modularisation.");
+                            isubdft.clear();
+                        } else {
+                            isubdft = {child->id()};
+                        }
+                    }
+                    if(isubdft.empty()) {
+                        subtreeBEs.clear();
+                        break;
+                    } else {
+                        // Add new independent subtree
+                        std::vector<size_t> subtree;
+                        for (size_t id : isubdft) {
+                            if (dft.isBasicElement(id)) {
+                                subtree.push_back(id);
+                            }
+                        }
+                        subtreeBEs.push_back(subtree);
+                    }
+                }
+            }
+            if (subtreeBEs.empty()) {
+                // Consider complete tree
+                std::vector<size_t> subtree;
+                for (size_t id = 0; id < dft.nrElements(); ++id) {
+                    if (dft.isBasicElement(id)) {
+                        subtree.push_back(id);
+                    }
+                }
+                subtreeBEs.push_back(subtree);
+            }
+
+            for (auto tree : subtreeBEs) {
+                std::stringstream stream;
+                stream << "Subtree: ";
+                for (size_t id : tree) {
+                    stream << id << ", ";
+                }
+                STORM_LOG_TRACE(stream.str());
+            }
         }
 
         template<typename ValueType, typename StateType>
@@ -476,18 +531,42 @@ namespace storm {
         template<typename ValueType, typename StateType>
         ValueType ExplicitDFTModelBuilderApprox<ValueType, StateType>::getUpperBound(DFTStatePointer const& state) const {
             // Get the upper bound by considering the failure of all BEs
-            // The used formula for the rate is 1/( 1/a + 1/b + 1/c + ... - 1/(a+b) - 1/(a+c) - ... + 1/(a+b+c) + ... - ...)
-            ValueType upperBound = storm::utility::zero<ValueType>();
+            ValueType upperBound = storm::utility::one<ValueType>();
+            ValueType rateSum = storm::utility::zero<ValueType>();
+
+            // Compute for each independent subtree
+            for (std::vector<size_t> const& subtree : subtreeBEs) {
+                // Get all possible rates
+                std::vector<ValueType> rates;
+                for (size_t id : subtree) {
+                    if (state->isOperational(id)) {
+                        // Get BE rate
+                        rates.push_back(state->getBERate(id, true));
+                        rateSum += rates.back();
+                    }
+                }
 
-            // Get all possible rates
-            std::vector<ValueType> rates(state->nrFailableBEs() + state->nrNotFailableBEs());
-            for (size_t index = 0; index < state->nrFailableBEs(); ++index) {
-                rates[index] = state->getFailableBERate(index);
+                // We move backwards and start with swapping the last element to itself
+                // Then we do not need to swap back
+                for (auto it = rates.rbegin(); it != rates.rend(); ++it) {
+                    // Compute AND MTTF of subtree without current rate and scale with current rate
+                    std::iter_swap(it, rates.end() - 1);
+                    upperBound += rates.back() * computeMTTFAnd(rates, rates.size() - 1);
+                }
             }
-            for (size_t index = 0; index < state->nrNotFailableBEs(); ++index) {
-                rates[index + state->nrFailableBEs()] = state->getNotFailableBERate(index);
+
+            STORM_LOG_TRACE("Upper bound is " << (rateSum / upperBound) << " for state " << state->getId());
+            STORM_LOG_ASSERT(!storm::utility::isZero(upperBound), "Upper bound is 0");
+            STORM_LOG_ASSERT(!storm::utility::isZero(rateSum), "State is absorbing");
+            return rateSum / upperBound;
+        }
+
+        template<typename ValueType, typename StateType>
+        ValueType ExplicitDFTModelBuilderApprox<ValueType, StateType>::computeMTTFAnd(std::vector<ValueType> rates, size_t size) const {
+            ValueType result = storm::utility::zero<ValueType>();
+            if (size == 0) {
+                return result;
             }
-            STORM_LOG_ASSERT(rates.size() > 0, "State is absorbing");
 
             // TODO Matthias: works only for <64 BEs!
             // WARNING: this code produces wrong results for more than 32 BEs
@@ -508,30 +587,30 @@ namespace storm {
                     sum += storm::utility::one<ValueType>() / permResult;
                 } while(permutation < (static_cast<size_t>(1) << rates.size()) && permutation != 0);
                 if (i % 2 == 0) {
-                    upperBound -= sum;
+                    result -= sum;
                 } else {
-                    upperBound += sum;
+                    result += sum;
                 }
             }*/
 
-            // Compute upper bound with permutations of size <= 3
-            for (size_t i1 = 0; i1 < rates.size(); ++i1) {
+            // Compute result with permutations of size <= 3
+            for (size_t i1 = 0; i1 < size; ++i1) {
                 // + 1/a
                 ValueType sum = rates[i1];
-                upperBound += storm::utility::one<ValueType>() / sum;
+                result += storm::utility::one<ValueType>() / sum;
                 for (size_t i2 = 0; i2 < i1; ++i2) {
                     // - 1/(a+b)
                     ValueType sum2 = sum + rates[i2];
-                    upperBound -= storm::utility::one<ValueType>() / sum2;
+                    result -= storm::utility::one<ValueType>() / sum2;
                     for (size_t i3 = 0; i3 < i2; ++i3) {
                         // + 1/(a+b+c)
-                        upperBound += storm::utility::one<ValueType>() / (sum2 + rates[i3]);
+                        result += storm::utility::one<ValueType>() / (sum2 + rates[i3]);
                     }
                 }
             }
 
-            STORM_LOG_ASSERT(!storm::utility::isZero(upperBound), "UpperBound is 0");
-            return storm::utility::one<ValueType>() / upperBound;
+            STORM_LOG_ASSERT(!storm::utility::isZero(result), "UpperBound is 0");
+            return result;
         }
 
         template<typename ValueType, typename StateType>

src/builder/ExplicitDFTModelBuilderApprox.h

@@ -238,6 +238,16 @@ namespace storm {
              */
             ValueType getUpperBound(DFTStatePointer const& state) const;
 
+            /*!
+             * Compute the MTTF of an AND gate via a closed formula.
+             * The used formula is 1/( 1/a + 1/b + 1/c + ... - 1/(a+b) - 1/(a+c) - ... + 1/(a+b+c) + ... - ...)
+             *
+             * @param rates List of rates of children of AND.
+             * @param size  Only indices < size are considered in the vector.
+             * @return MTTF.
+             */
+            ValueType computeMTTFAnd(std::vector<ValueType> rates, size_t size) const;
+
             /*!
              * Compares the priority of two states.
              *
@@ -314,6 +324,9 @@ namespace storm {
             // Notice that we need an ordered map here to easily iterate in increasing order over state ids.
             // TODO remove again
             std::map<StateType, std::pair<DFTStatePointer, ExplorationHeuristicPointer>> skippedStates;
+
+            // List of independent subtrees and the BEs contained in them.
+            std::vector<std::vector<size_t>> subtreeBEs;
         };
 
     }

src/storage/BucketPriorityQueue.cpp

@@ -41,7 +41,6 @@ namespace storm {
                 return immediateBucket.back();
             }
             STORM_LOG_ASSERT(!empty(), "BucketPriorityQueue is empty");
-            STORM_LOG_ASSERT(nrUnsortedItems == 0, "First bucket is not sorted");
             return buckets[currentBucket].front();
         }
 
@@ -126,7 +125,6 @@ namespace storm {
                 return;
             }
             STORM_LOG_ASSERT(!empty(), "BucketPriorityQueue is empty");
-            STORM_LOG_ASSERT(nrUnsortedItems == 0, "First bucket is not sorted");
             std::pop_heap(buckets[currentBucket].begin(), buckets[currentBucket].end(), compare);
             buckets[currentBucket].pop_back();
             if (buckets[currentBucket].empty()) {

src/storage/dft/DFTState.cpp

@@ -241,21 +241,26 @@ namespace storm {
         }
 
         template<typename ValueType>
-        ValueType DFTState<ValueType>::getBERate(size_t id, bool considerPassive) const {
+        ValueType DFTState<ValueType>::getBERate(size_t id, bool avoidCold) const {
             STORM_LOG_ASSERT(mDft.isBasicElement(id), "Element is no BE.");
-            if (considerPassive && mDft.hasRepresentant(id) && !isActive(mDft.getRepresentant(id))) {
-                STORM_LOG_ASSERT(!storm::utility::isZero(mDft.getBasicElement(id)->passiveFailureRate()), "Failure rate of BE " << mDft.getBasicElement(id)->id() << " is 0 in state " << mDft.getStateString(mStatus, mStateGenerationInfo, mId));
-                return mDft.getBasicElement(id)->passiveFailureRate();
-            } else {
-                STORM_LOG_ASSERT(!storm::utility::isZero(mDft.getBasicElement(id)->activeFailureRate()), "Failure rate of BE " << mDft.getBasicElement(id)->id() << " is 0 in state " << mDft.getStateString(mStatus, mStateGenerationInfo, mId));
-                return mDft.getBasicElement(id)->activeFailureRate();
+
+            if (mDft.hasRepresentant(id) && !isActive(mDft.getRepresentant(id))) {
+                if (avoidCold && mDft.getBasicElement(id)->isColdBasicElement()) {
+                    // Return active failure rate instead of 0.
+                    return mDft.getBasicElement(id)->activeFailureRate();
+                } else {
+                    // Return passive failure rate
+                    return mDft.getBasicElement(id)->passiveFailureRate();
+                }
             }
+            // Return active failure rate
+            return mDft.getBasicElement(id)->activeFailureRate();
         }
 
         template<typename ValueType>
         ValueType DFTState<ValueType>::getFailableBERate(size_t index) const {
             STORM_LOG_ASSERT(index < nrFailableBEs(), "Index invalid.");
-            return getBERate(mCurrentlyFailableBE[index], true);
+            return getBERate(mCurrentlyFailableBE[index], false);
         }
 
         template<typename ValueType>
@@ -264,7 +269,7 @@ namespace storm {
             STORM_LOG_ASSERT(isPseudoState() || storm::utility::isZero<ValueType>(mDft.getBasicElement(mCurrentlyNotFailableBE[index])->activeFailureRate()) ||
                              (mDft.hasRepresentant(mCurrentlyNotFailableBE[index]) && !isActive(mDft.getRepresentant(mCurrentlyNotFailableBE[index]))), "BE " << mCurrentlyNotFailableBE[index] << " can fail in state: " << mDft.getStateString(mStatus, mStateGenerationInfo, mId));
             // Use active failure rate as passive failure rate is 0.
-            return getBERate(mCurrentlyNotFailableBE[index], false);
+            return getBERate(mCurrentlyNotFailableBE[index], true);
         }
 
         template<typename ValueType>

src/storage/dft/DFTState.h

@@ -208,6 +208,17 @@ namespace storm {
              */
             ValueType getNotFailableBERate(size_t index) const;
 
+            /**
+             * Get the failure rate of the given BE.
+             *
+             * @param id        Id of BE.
+             * @param avoidCold Flag indicating if a cold passive failure rate should be avoided by giving
+             *                  the active failure rate instead.
+             *
+             * @return Failure rate of the BE.
+             */
+            ValueType getBERate(size_t id, bool avoidCold) const;
+
             /** Get number of currently failable dependencies.
              *
              * @return Number of failable dependencies.
@@ -302,17 +313,6 @@ namespace storm {
         private:
             void propagateActivation(size_t representativeId);
 
-            /**
-             * Get the failure rate of the given BE.
-             *
-             * @param id              Id of BE.
-             * @param considerPassive Flag indicating if the passive failure rate should be returned if
-             *                        the BE currently is passive.
-             *
-             * @return Failure rate of the BE.
-             */
-            ValueType getBERate(size_t id, bool considerPassive) const;
-
         };
 
     }