ValueIterationOperator.h

Go to the documentation of this file.

#pragma once
#include <functional>
#include <optional>
#include <utility>
#include <vector>
 
#include <boost/range/adaptor/reversed.hpp>
#include <boost/range/irange.hpp>
 
#include "storm/solver/helper/ValueIterationOperatorForward.h"
#include "storm/storage/sparse/StateType.h"
#include "storm/utility/macros.h"
#include "storm/utility/vector.h"  // TODO
 
namespace storm {
class Environment;
 
namespace storage {
template<typename T>
class SparseMatrix;
}
 
namespace solver::helper {
 
template<typename ValueType, bool TrivialRowGrouping, typename SolutionType>
class ValueIterationOperator {
   public:
    using IndexType = storm::storage::sparse::state_type;
 
    template<bool Backward = true>
    void setMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
 
    void setMatrixForwards(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
 
    void setMatrixBackwards(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
 
    template<typename OperandType, typename OffsetType, typename BackendType>
    bool apply(OperandType const& operandIn, OperandType& operandOut, OffsetType const& offsets, BackendType& backend) const {
        // We assume in the normal apply case that we can just maximize.
        return applyRobust<OptimizationDirection::Maximize>(operandIn, operandOut, offsets, backend);
    }
 
    template<OptimizationDirection RobustDir, typename OperandType, typename OffsetType, typename BackendType>
    bool applyRobust(OperandType const& operandIn, OperandType& operandOut, OffsetType const& offsets, BackendType& backend) const {
        if (hasSkippedRows) {
            if (backwards) {
                return apply<OperandType, OffsetType, BackendType, true, true, RobustDir>(operandOut, operandIn, offsets, backend);
            } else {
                return apply<OperandType, OffsetType, BackendType, false, true, RobustDir>(operandOut, operandIn, offsets, backend);
            }
        } else {
            if (backwards) {
                return apply<OperandType, OffsetType, BackendType, true, false, RobustDir>(operandOut, operandIn, offsets, backend);
            } else {
                return apply<OperandType, OffsetType, BackendType, false, false, RobustDir>(operandOut, operandIn, offsets, backend);
            }
        }
    }
 
    template<typename OperandType, typename OffsetType, typename BackendType>
    bool applyInPlace(OperandType& operand, OffsetType const& offsets, BackendType& backend) const {
        return apply(operand, operand, offsets, backend);
    }
 
    template<OptimizationDirection RobustDir, typename OperandType, typename OffsetType, typename BackendType>
    bool applyInPlaceRobust(OperandType& operand, OffsetType const& offsets, BackendType& backend) const {
        return applyRobust<RobustDir>(operand, operand, offsets, backend);
    }
 
    void setIgnoredRows(bool useLocalRowIndices, std::function<bool(IndexType, IndexType)> const& ignore);
 
    void unsetIgnoredRows();
 
    std::vector<IndexType> const& getRowGroupIndices() const;
 
    std::vector<SolutionType>& allocateAuxiliaryVector(uint64_t size, std::optional<SolutionType> const& initialValue = {});
 
    void freeAuxiliaryVector();
 
   private:
    template<typename OperandType, typename OffsetType, typename BackendType, bool Backward, bool SkipIgnoredRows, OptimizationDirection RobustDirection>
    bool apply(OperandType& operandOut, OperandType const& operandIn, OffsetType const& offsets, BackendType& backend) const {
        STORM_LOG_ASSERT(getSize(operandIn) == getSize(operandOut), "Input and Output Operands have different sizes.");
        auto const operandSize = getSize(operandIn);
        STORM_LOG_ASSERT(TrivialRowGrouping || rowGroupIndices->size() == operandSize + 1, "Dimension mismatch");
        backend.startNewIteration();
        auto matrixValueIt = matrixValues.cbegin();
        auto matrixColumnIt = matrixColumns.cbegin();
        for (auto groupIndex : indexRange<Backward>(0, operandSize)) {
            STORM_LOG_ASSERT(matrixColumnIt != matrixColumns.end(), "VI Operator in invalid state.");
            STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
            //            STORM_LOG_ASSERT(matrixValueIt != matrixValues.end(), "VI Operator in invalid state.");
            if constexpr (TrivialRowGrouping) {
                backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, groupIndex), groupIndex, groupIndex);
            } else {
                IndexType rowIndex = (*rowGroupIndices)[groupIndex];
                if constexpr (SkipIgnoredRows) {
                    rowIndex += skipMultipleIgnoredRows(matrixColumnIt, matrixValueIt);
                }
                backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
                while (*matrixColumnIt < StartOfRowGroupIndicator) {
                    ++rowIndex;
                    if (!SkipIgnoredRows || !skipIgnoredRow(matrixColumnIt, matrixValueIt)) {
                        backend.nextRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
                    }
                }
            }
            if constexpr (isPair<OperandType>::value) {
                backend.applyUpdate(operandOut.first[groupIndex], operandOut.second[groupIndex], groupIndex);
            } else {
                backend.applyUpdate(operandOut[groupIndex], groupIndex);
            }
            if (backend.abort()) {
                return backend.converged();
            }
        }
        STORM_LOG_ASSERT(matrixColumnIt + 1 == matrixColumns.cend(), "Unexpected position of matrix column iterator.");
        STORM_LOG_ASSERT(matrixValueIt == matrixValues.cend(), "Unexpected position of matrix column iterator.");
        backend.endOfIteration();
        return backend.converged();
    }
 
    // Auxiliary methods to deal with various OperandTypes and OffsetTypes
 
    template<typename OpT, typename OffT>
    OpT initializeRowRes(std::vector<OpT> const&, std::vector<OffT> const& offsets, uint64_t offsetIndex) const {
        return offsets[offsetIndex];
    }
 
    template<typename OpT1, typename OpT2, typename OffT>
    std::pair<OpT1, OpT2> initializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::vector<OffT> const& offsets,
                                           uint64_t offsetIndex) const {
        return {offsets[offsetIndex], offsets[offsetIndex]};
    }
 
    template<typename OpT1, typename OpT2, typename OffT1, typename OffT2>
    std::pair<OpT1, OpT2> initializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::pair<std::vector<OffT1> const*, OffT2> const& offsets,
                                           uint64_t offsetIndex) const {
        return {(*offsets.first)[offsetIndex], offsets.second};
    }
 
    template<OptimizationDirection RobustDirection, typename OpT, typename OffT>
    OpT robustInitializeRowRes(std::vector<OpT> const&, std::vector<OffT> const& offsets, uint64_t offsetIndex) const {
        return offsets[offsetIndex].upper();
    }
 
    template<OptimizationDirection RobustDirection, typename OpT1, typename OpT2, typename OffT>
    std::pair<OpT1, OpT2> robustInitializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::vector<OffT> const& offsets,
                                                 uint64_t offsetIndex) const {
        STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
 
        return {offsets[offsetIndex].upper(), offsets[offsetIndex].upper()};
    }
 
    template<OptimizationDirection RobustDirection, typename OpT1, typename OpT2, typename OffT1, typename OffT2>
    std::pair<OpT1, OpT2> robustInitializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&,
                                                 std::pair<std::vector<OffT1> const*, OffT2> const& offsets, uint64_t offsetIndex) const {
        STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
 
        return {(*offsets.first)[offsetIndex], offsets.second};
    }
 
    template<OptimizationDirection RobustDirection, typename OperandType, typename OffsetType>
    auto applyRow(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
                  OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
        if constexpr (std::is_same_v<ValueType, storm::Interval>) {
            return applyRowRobust<RobustDirection>(matrixColumnIt, matrixValueIt, operand, offsets, offsetIndex);
        } else {
            return applyRowStandard(matrixColumnIt, matrixValueIt, operand, offsets, offsetIndex);
        }
    }
 
    template<typename OperandType, typename OffsetType>
    auto applyRowStandard(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
                          OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
        STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
        auto result{initializeRowRes(operand, offsets, offsetIndex)};
        for (++matrixColumnIt; *matrixColumnIt < StartOfRowIndicator; ++matrixColumnIt, ++matrixValueIt) {
            if constexpr (isPair<OperandType>::value) {
                result.first += operand.first[*matrixColumnIt] * (*matrixValueIt);
                result.second += operand.second[*matrixColumnIt] * (*matrixValueIt);
            } else {
                result += operand[*matrixColumnIt] * (*matrixValueIt);
            }
        }
        return result;
    }
 
    // Aux function for applyRowRobust
    template<OptimizationDirection RobustDirection>
    struct AuxCompare {
        bool operator()(const std::pair<SolutionType, SolutionType>& a, const std::pair<SolutionType, SolutionType>& b) const {
            if constexpr (RobustDirection == OptimizationDirection::Maximize) {
                return a.first > b.first;
            } else {
                return a.first < b.first;
            }
        }
    };
 
    template<OptimizationDirection RobustDirection, typename OperandType, typename OffsetType>
    auto applyRowRobust(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
                        OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
        STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
        auto result{robustInitializeRowRes<RobustDirection>(operand, offsets, offsetIndex)};
        AuxCompare<RobustDirection> compare;
        applyCache.robustOrder.clear();
 
        SolutionType remainingValue{storm::utility::one<SolutionType>()};
        for (++matrixColumnIt; *matrixColumnIt < StartOfRowIndicator; ++matrixColumnIt, ++matrixValueIt) {
            auto const lower = matrixValueIt->lower();
            if constexpr (isPair<OperandType>::value) {
                STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
                // Notice the unclear semantics here in terms of how to order things.
            } else {
                result += operand[*matrixColumnIt] * lower;
            }
            remainingValue -= lower;
            auto const diameter = matrixValueIt->upper() - lower;
            if (!storm::utility::isZero(diameter)) {
                applyCache.robustOrder.emplace_back(operand[*matrixColumnIt], diameter);
            }
        }
        if (storm::utility::isZero(remainingValue) || storm::utility::isOne(remainingValue)) {
            return result;
        }
 
        std::sort(applyCache.robustOrder.begin(), applyCache.robustOrder.end(), compare);
 
        for (auto const& pair : applyCache.robustOrder) {
            auto availableMass = std::min(pair.second, remainingValue);
            result += availableMass * pair.first;
            remainingValue -= availableMass;
            if (storm::utility::isZero(remainingValue)) {
                return result;
            }
        }
        STORM_LOG_ASSERT(storm::utility::isAlmostZero(remainingValue), "Remaining value should be zero (all prob mass taken) but is " << remainingValue);
        return result;
    }
 
    // Auxiliary helpers used for metaprogramming
    template<bool Backward>
    auto indexRange(IndexType start, IndexType end) const {
        if constexpr (Backward) {
            return boost::adaptors::reverse(boost::irange(start, end));
        } else {
            return boost::irange(start, end);
        }
    }
 
    template<typename T>
    uint64_t getSize(std::vector<T> const& vec) const {
        return vec.size();
    }
 
    template<typename T1, typename T2>
    uint64_t getSize(std::pair<T1, T2> const& pairOfVec) const {
        return pairOfVec.first.size();
    }
 
    template<typename>
    struct isPair : std::false_type {};
 
    template<typename T1, typename T2>
    struct isPair<std::pair<T1, T2>> : std::true_type {};
 
    template<bool Backward = true>
    void setIgnoredRows(bool useLocalRowIndices, std::function<bool(IndexType, IndexType)> const& ignore);
 
    void moveToEndOfRow(std::vector<IndexType>::iterator& matrixColumnIt) const;
 
    bool skipIgnoredRow(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt) const;
 
    uint64_t skipMultipleIgnoredRows(std::vector<IndexType>::const_iterator& matrixColumnIt,
                                     typename std::vector<ValueType>::const_iterator& matrixValueIt) const;
 
    std::vector<ValueType> matrixValues;
 
    std::vector<IndexType> matrixColumns;
 
    std::vector<IndexType> const* rowGroupIndices;
 
    bool backwards{true};
 
    bool hasSkippedRows{false};
 
    std::vector<SolutionType> auxiliaryVector;
 
    bool auxiliaryVectorUsedExternally{false};
 
    // Due to a GCC bug we have to add this dummy template type here
    // https://stackoverflow.com/questions/49707184/explicit-specialization-in-non-namespace-scope-does-not-compile-in-gcc
    template<typename ApplyValueType, typename Dummy>
    struct ApplyCache {};
 
    template<typename Dummy>
    struct ApplyCache<storm::Interval, Dummy> {
        mutable std::vector<std::pair<SolutionType, SolutionType>> robustOrder;
    };
 
    ApplyCache<ValueType, int> applyCache;
 
    IndexType const StartOfRowIndicator = 1ull << 63;  // 10000..0
 
    IndexType const StartOfRowGroupIndicator = StartOfRowIndicator + (1ull << 62);  // 11000..0
 
    IndexType const SkipNumEntriesMask = ~StartOfRowGroupIndicator;  // 00111..1
};
 
}  // namespace solver::helper
}  // namespace storm