vector.h

Go to the documentation of this file.

#ifndef STORM_UTILITY_VECTOR_H_
#define STORM_UTILITY_VECTOR_H_
 
#include <algorithm>
#include <functional>
#include <iosfwd>
#include <numeric>
#include "storm/adapters/IntelTbbAdapter.h"
#include "storm/adapters/RationalNumberAdapter.h"
 
#include <boost/optional.hpp>
 
#include "storm/solver/OptimizationDirection.h"
#include "storm/storage/BitVector.h"
#include "storm/utility/constants.h"
#include "storm/utility/macros.h"
 
#include "storm/exceptions/NotImplementedException.h"
 
namespace storm {
namespace utility {
namespace vector {
 
template<typename ValueType>
struct VectorHash {
    size_t operator()(std::vector<ValueType> const& vec) const {
        std::hash<ValueType> hasher;
        std::size_t seed = 0;
        for (ValueType const& element : vec) {
            seed ^= hasher(element) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
        }
        return seed;
    }
};
 
template<class T>
std::size_t findOrInsert(std::vector<T>& vector, T&& element) {
    std::size_t position = std::find(vector.begin(), vector.end(), element) - vector.begin();
    if (position == vector.size()) {
        vector.emplace_back(std::move(element));
    }
    return position;
}
 
template<typename T>
void setAllValues(std::vector<T>& vec, storm::storage::BitVector const& positions, T const& positiveValue = storm::utility::one<T>(),
                  T const& negativeValue = storm::utility::zero<T>()) {
    if (positions.getNumberOfSetBits() * 2 > positions.size()) {
        vec.resize(positions.size(), positiveValue);
        uint64_t index = positions.getNextUnsetIndex(0);
        while (index < positions.size()) {
            vec[index] = negativeValue;
            index = positions.getNextUnsetIndex(index + 1);
        }
    } else {
        vec.resize(positions.size(), negativeValue);
        for (uint64_t index : positions) {
            vec[index] = positiveValue;
        }
    }
}
 
template<class T>
void setVectorValues(std::vector<T>& vector, storm::storage::BitVector const& positions, std::vector<T> const& values) {
    STORM_LOG_ASSERT(positions.size() <= vector.size(), "We cannot set positions that have not been initialized");
    STORM_LOG_ASSERT(positions.getNumberOfSetBits() <= values.size(), "The number of selected positions (" << positions.getNumberOfSetBits()
                                                                                                           << ") exceeds the size of the input vector ("
                                                                                                           << values.size() << ").");
    uint_fast64_t oldPosition = 0;
    for (auto position : positions) {
        vector[position] = values[oldPosition++];
    }
}
 
template<class T>
void setVectorValues(std::vector<T>& vector, storm::storage::BitVector const& positions, T value) {
    STORM_LOG_ASSERT(positions.size() <= vector.size(), "We cannot set positions that have not been initialized");
    for (auto position : positions) {
        vector[position] = value;
    }
}
 
template<typename T>
void setNonzeroIndices(std::vector<T> const& vec, storm::storage::BitVector& bv) {
    STORM_LOG_ASSERT(bv.size() == vec.size(), "Bitvector size should match vector size");
    for (uint64_t i = 0; i < vec.size(); ++i) {
        if (!storm::utility::isZero(vec[i])) {
            bv.set(i, true);
        }
    }
}
 
template<class OutputIterator, class Size, class Assignable>
void iota_n(OutputIterator first, Size n, Assignable value) {
    std::generate_n(first, n, [&value]() { return value++; });
}
 
template<typename T>
inline std::vector<T> buildVectorForRange(T min, T max) {
    STORM_LOG_ASSERT(min <= max, "Invalid range.");
    T diff = max - min;
    std::vector<T> v;
    v.reserve(diff);
    iota_n(std::back_inserter(v), diff, min);
    return v;
}
 
template<typename T>
std::vector<uint_fast64_t> getSortedIndices(std::vector<T> const& v) {
    std::vector<uint_fast64_t> res = buildVectorForRange<uint_fast64_t>(0, v.size());
    std::sort(res.begin(), res.end(), [&v](uint_fast64_t index1, uint_fast64_t index2) { return v[index1] > v[index2]; });
    return res;
}
 
template<typename T>
bool isUnique(std::vector<T> const& v) {
    if (v.size() < 2) {
        return true;
    }
    auto sortedIndices = getSortedIndices(v);
    auto indexIt = sortedIndices.begin();
    T const* previous = &v[*indexIt];
    for (++indexIt; indexIt != sortedIndices.end(); ++indexIt) {
        T const& current = v[*indexIt];
        if (current == *previous) {
            return false;
        }
        previous = &current;
    }
    return true;
}
 
template<typename T, typename Comparator>
bool compareElementWise(std::vector<T> const& left, std::vector<T> const& right, Comparator comp = std::less<T>()) {
    STORM_LOG_ASSERT(left.size() == right.size(), "Expected that vectors for comparison have equal size");
    return std::equal(left.begin(), left.end(), right.begin(), comp);
}
 
template<class T>
void selectVectorValues(std::vector<T>& vector, storm::storage::BitVector const& positions, std::vector<T> const& values) {
    STORM_LOG_ASSERT(positions.getNumberOfSetBits() <= vector.size(), "The number of selected positions (" << positions.getNumberOfSetBits()
                                                                                                           << ") exceeds the size of the target vector ("
                                                                                                           << vector.size() << ").");
    STORM_LOG_ASSERT(positions.size() == values.size(),
                     "Size mismatch of the positions vector (" << positions.size() << ") and the values vector (" << values.size() << ").");
    auto targetIt = vector.begin();
    for (auto position : positions) {
        *targetIt = values[position];
        ++targetIt;
    }
}
 
template<class T>
void selectVectorValues(std::vector<T>& vector, storm::storage::BitVector const& positions, std::vector<uint_fast64_t> const& rowGrouping,
                        std::vector<T> const& values) {
    auto targetIt = vector.begin();
    for (auto position : positions) {
        for (uint_fast64_t i = rowGrouping[position]; i < rowGrouping[position + 1]; ++i, ++targetIt) {
            *targetIt = values[i];
        }
    }
}
 
template<class T>
void selectVectorValues(std::vector<T>& vector, std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGrouping,
                        std::vector<T> const& values) {
    auto targetIt = vector.begin();
    for (uint_fast64_t i = 0; i < vector.size(); ++i, ++targetIt) {
        *targetIt = values[rowGrouping[i] + rowGroupToRowIndexMapping[i]];
    }
}
 
template<class T>
void selectVectorValues(std::vector<T>& vector, std::vector<uint_fast64_t> const& indexSequence, std::vector<T> const& values) {
    STORM_LOG_ASSERT(indexSequence.size() <= vector.size(),
                     "The number of selected positions (" << indexSequence.size() << ") exceeds the size of the target vector (" << vector.size() << ").");
 
    for (uint_fast64_t vectorIndex = 0; vectorIndex < vector.size(); ++vectorIndex) {
        vector[vectorIndex] = values[indexSequence[vectorIndex]];
    }
}
 
template<class T>
void selectVectorValuesRepeatedly(std::vector<T>& vector, storm::storage::BitVector const& positions, std::vector<uint_fast64_t> const& rowGrouping,
                                  std::vector<T> const& values) {
    auto targetIt = vector.begin();
    for (auto position : positions) {
        for (uint_fast64_t i = rowGrouping[position]; i < rowGrouping[position + 1]; ++i, ++targetIt) {
            *targetIt = values[position];
        }
    }
}
 
template<class T>
void subtractFromConstantOneVector(std::vector<T>& vector) {
    for (auto& element : vector) {
        element = storm::utility::one<T>() - element;
    }
}
 
template<class T>
void addFilteredVectorGroupsToGroupedVector(std::vector<T>& target, std::vector<T> const& source, storm::storage::BitVector const& filter,
                                            std::vector<uint_fast64_t> const& rowGroupIndices) {
    auto targetIt = target.begin();
    for (auto group : filter) {
        auto it = source.cbegin() + rowGroupIndices[group];
        auto ite = source.cbegin() + rowGroupIndices[group + 1];
        for (; it != ite; ++targetIt, ++it) {
            *targetIt += *it;
        }
    }
}
 
template<class T>
void addVectorToGroupedVector(std::vector<T>& target, std::vector<T> const& source, std::vector<uint_fast64_t> const& rowGroupIndices) {
    auto targetIt = target.begin();
    auto sourceIt = source.cbegin();
    auto sourceIte = source.cend();
    auto rowGroupIt = rowGroupIndices.cbegin();
 
    for (; sourceIt != sourceIte; ++sourceIt) {
        uint_fast64_t current = *rowGroupIt;
        ++rowGroupIt;
        uint_fast64_t next = *rowGroupIt;
        for (; current < next; ++targetIt, ++current) {
            *targetIt += *sourceIt;
        }
    }
}
 
template<class T>
void addFilteredVectorToGroupedVector(std::vector<T>& target, std::vector<T> const& source, storm::storage::BitVector const& filter,
                                      std::vector<uint_fast64_t> const& rowGroupIndices) {
    auto targetIt = target.begin();
    for (auto group : filter) {
        uint_fast64_t current = rowGroupIndices[group];
        uint_fast64_t next = rowGroupIndices[group + 1];
        for (; current < next; ++current, ++targetIt) {
            *targetIt += source[group];
        }
    }
}
 
template<class InValueType1, class InValueType2, class OutValueType, class Operation>
void applyPointwiseTernary(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target,
                           Operation f = Operation()) {
    auto firstIt = firstOperand.begin();
    auto firstIte = firstOperand.end();
    auto secondIt = secondOperand.begin();
    auto targetIt = target.begin();
    while (firstIt != firstIte) {
        *targetIt = f(*firstIt, *secondIt, *targetIt);
        ++targetIt;
        ++firstIt;
        ++secondIt;
    }
}
 
#ifdef STORM_HAVE_INTELTBB
template<class InValueType1, class InValueType2, class OutValueType, class Operation>
void applyPointwiseTernaryParallel(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand,
                                   std::vector<OutValueType>& target, Operation f = Operation()) {
    tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()), [&](tbb::blocked_range<uint_fast64_t> const& range) {
        auto firstIt = firstOperand.begin() + range.begin();
        auto firstIte = firstOperand.begin() + range.end();
        auto secondIt = secondOperand.begin() + range.begin();
        auto targetIt = target.begin() + range.begin();
        while (firstIt != firstIte) {
            *targetIt = f(*firstIt, *secondIt, *targetIt);
            ++targetIt;
            ++firstIt;
            ++secondIt;
        }
    });
}
#endif
 
template<class InValueType1, class InValueType2, class OutValueType, class Operation>
void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target,
                    Operation f = Operation()) {
    std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), f);
}
 
#ifdef STORM_HAVE_INTELTBB
template<class InValueType1, class InValueType2, class OutValueType, class Operation>
void applyPointwiseParallel(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target,
                            Operation f = Operation()) {
    tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()), [&](tbb::blocked_range<uint_fast64_t> const& range) {
        std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(),
                       target.begin() + range.begin(), f);
    });
}
#endif
 
template<class InValueType, class OutValueType, class Operation>
void applyPointwise(std::vector<InValueType> const& operand, std::vector<OutValueType>& target, Operation f = Operation()) {
    std::transform(operand.begin(), operand.end(), target.begin(), f);
}
 
#ifdef STORM_HAVE_INTELTBB
template<class InValueType, class OutValueType, class Operation>
void applyPointwiseParallel(std::vector<InValueType> const& operand, std::vector<OutValueType>& target, Operation f = Operation()) {
    tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()), [&](tbb::blocked_range<uint_fast64_t> const& range) {
        std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), f);
    });
}
#endif
 
template<class InValueType1, class InValueType2, class OutValueType>
void addVectors(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
    applyPointwise<InValueType1, InValueType2, OutValueType, std::plus<>>(firstOperand, secondOperand, target);
}
 
template<class InValueType1, class InValueType2, class OutValueType>
void subtractVectors(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
    applyPointwise<InValueType1, InValueType2, OutValueType, std::minus<>>(firstOperand, secondOperand, target);
}
 
template<class InValueType1, class InValueType2, class OutValueType>
void multiplyVectorsPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand,
                              std::vector<OutValueType>& target) {
    applyPointwise<InValueType1, InValueType2, OutValueType, std::multiplies<>>(firstOperand, secondOperand, target);
}
 
template<class InValueType1, class InValueType2, class OutValueType>
void divideVectorsPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target) {
    applyPointwise<InValueType1, InValueType2, OutValueType, std::divides<>>(firstOperand, secondOperand, target);
}
 
template<class ValueType1, class ValueType2>
void scaleVectorInPlace(std::vector<ValueType1>& target, ValueType2 const& factor) {
    applyPointwise<ValueType1, ValueType1>(target, target, [&](ValueType1 const& argument) -> ValueType1 { return argument * factor; });
}
 
template<class InValueType1, class InValueType2, class InValueType3>
void addScaledVector(std::vector<InValueType1>& firstOperand, std::vector<InValueType2> const& secondOperand, InValueType3 const& factor) {
    applyPointwise<InValueType1, InValueType2, InValueType1>(
        firstOperand, secondOperand, firstOperand, [&](InValueType1 const& val1, InValueType2 const& val2) -> InValueType1 { return val1 + (factor * val2); });
}
 
template<class T>
T dotProduct(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand) {
    return std::inner_product(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), storm::utility::zero<T>());
}
 
template<class T>
storm::storage::BitVector filter(std::vector<T> const& values, std::function<bool(T const& value)> const& function) {
    storm::storage::BitVector result(values.size(), false);
 
    uint_fast64_t currentIndex = 0;
    for (auto const& value : values) {
        if (function(value)) {
            result.set(currentIndex, true);
        }
        ++currentIndex;
    }
 
    return result;
}
 
template<class T>
storm::storage::BitVector filterGreaterZero(std::vector<T> const& values) {
    return filter<T>(values, [](T const& value) -> bool { return value > storm::utility::zero<T>(); });
}
 
template<class T>
storm::storage::BitVector filterZero(std::vector<T> const& values) {
    return filter<T>(values, storm::utility::isZero<T>);
}
 
template<class T>
storm::storage::BitVector filterOne(std::vector<T> const& values) {
    return filter<T>(values, storm::utility::isOne<T>);
}
 
template<class T>
storm::storage::BitVector filterInfinity(std::vector<T> const& values) {
    return filter<T>(values, storm::utility::isInfinity<T>);
}
 
template<typename VT>
VT sum_if(std::vector<VT> const& values, storm::storage::BitVector const& filter) {
    STORM_LOG_ASSERT(values.size() == filter.size(), "Vector sizes mismatch.");
    VT sum = storm::utility::zero<VT>();
    for (auto pos : filter) {
        sum += values[pos];
    }
    return sum;
}
 
template<typename VT>
VT max_if(std::vector<VT> const& values, storm::storage::BitVector const& filter) {
    STORM_LOG_ASSERT(values.size() == filter.size(), "Vector sizes mismatch.");
    STORM_LOG_ASSERT(!filter.empty(), "Empty selection.");
 
    auto it = filter.begin();
    auto ite = filter.end();
 
    VT current = values[*it];
    ++it;
 
    for (; it != ite; ++it) {
        current = std::max(values[*it], current);
    }
    return current;
}
 
template<typename VT>
VT min_if(std::vector<VT> const& values, storm::storage::BitVector const& filter) {
    STORM_LOG_ASSERT(values.size() == filter.size(), "Vector sizes mismatch.");
    STORM_LOG_ASSERT(!filter.empty(), "Empty selection.");
 
    auto it = filter.begin();
    auto ite = filter.end();
 
    VT current = values[*it];
    ++it;
 
    for (; it != ite; ++it) {
        current = std::min(values[*it], current);
    }
    return current;
}
 
#ifdef STORM_HAVE_INTELTBB
template<class T, class Filter>
class TbbReduceVectorFunctor {
   public:
    TbbReduceVectorFunctor(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                           std::vector<uint_fast64_t>* choices, Filter const& f)
        : source(source), target(target), rowGrouping(rowGrouping), choices(choices), f(f) {
        // Intentionally left empty.
    }
 
    void operator()(tbb::blocked_range<uint64_t> const& range) const {
        uint_fast64_t startRow = range.begin();
        uint_fast64_t endRow = range.end();
 
        typename std::vector<T>::iterator targetIt = target.begin() + startRow;
        typename std::vector<T>::iterator targetIte = target.begin() + endRow;
        typename std::vector<uint_fast64_t>::const_iterator rowGroupingIt = rowGrouping.begin() + startRow;
        typename std::vector<T>::const_iterator sourceIt = source.begin() + *rowGroupingIt;
        typename std::vector<T>::const_iterator sourceIte;
        typename std::vector<uint_fast64_t>::iterator choiceIt;
        if (choices) {
            choiceIt = choices->begin() + startRow;
        }
 
        // Variables for correctly tracking choices (only update if new choice is strictly better).
        T oldSelectedChoiceValue;
        uint64_t selectedChoice;
 
        uint64_t currentRow = 0;
        for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) {
            // Only traverse elements if the row group is non-empty.
            if (*rowGroupingIt != *(rowGroupingIt + 1)) {
                *targetIt = *sourceIt;
 
                if (choices) {
                    selectedChoice = 0;
                    if (*choiceIt == 0) {
                        oldSelectedChoiceValue = *targetIt;
                    }
                }
 
                ++sourceIt;
                ++currentRow;
 
                for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++currentRow) {
                    if (choices && *choiceIt + *rowGroupingIt == currentRow) {
                        oldSelectedChoiceValue = *sourceIt;
                    }
 
                    if (f(*sourceIt, *targetIt)) {
                        *targetIt = *sourceIt;
                        if (choices) {
                            selectedChoice = std::distance(source.begin(), sourceIt) - *rowGroupingIt;
                        }
                    }
                }
 
                if (choices && f(*targetIt, oldSelectedChoiceValue)) {
                    *choiceIt = selectedChoice;
                }
            }
        }
    }
 
   private:
    std::vector<T> const& source;
    std::vector<T>& target;
    std::vector<uint_fast64_t> const& rowGrouping;
    std::vector<uint_fast64_t>* choices;
    Filter const& f;
};
#endif
 
template<class T, class Filter>
void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices) {
    Filter f;
    typename std::vector<T>::iterator targetIt = target.begin();
    typename std::vector<T>::iterator targetIte = target.end();
    typename std::vector<uint_fast64_t>::const_iterator rowGroupingIt = rowGrouping.begin();
    typename std::vector<T>::const_iterator sourceIt = source.begin();
    typename std::vector<T>::const_iterator sourceIte;
    typename std::vector<uint_fast64_t>::iterator choiceIt;
    if (choices) {
        choiceIt = choices->begin();
    }
 
    // Variables for correctly tracking choices (only update if new choice is strictly better).
    T oldSelectedChoiceValue;
    uint64_t selectedChoice;
 
    uint64_t currentRow = 0;
    for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) {
        // Only traverse elements if the row group is non-empty.
        if (*rowGroupingIt != *(rowGroupingIt + 1)) {
            *targetIt = *sourceIt;
 
            if (choices) {
                selectedChoice = 0;
                if (*choiceIt == 0) {
                    oldSelectedChoiceValue = *targetIt;
                }
            }
 
            ++sourceIt;
            ++currentRow;
 
            for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++currentRow) {
                if (choices && *rowGroupingIt + *choiceIt == currentRow) {
                    oldSelectedChoiceValue = *sourceIt;
                }
 
                if (f(*sourceIt, *targetIt)) {
                    *targetIt = *sourceIt;
                    if (choices) {
                        selectedChoice = std::distance(source.begin(), sourceIt) - *rowGroupingIt;
                    }
                }
            }
 
            if (choices && f(*targetIt, oldSelectedChoiceValue)) {
                *choiceIt = selectedChoice;
            }
        } else {
            *choiceIt = 0;
            *targetIt = storm::utility::zero<T>();
        }
    }
}
 
#ifdef STORM_HAVE_INTELTBB
template<class T, class Filter>
void reduceVectorParallel(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                          std::vector<uint_fast64_t>* choices) {
    tbb::parallel_for(tbb::blocked_range<uint64_t>(0, target.size()), TbbReduceVectorFunctor<T, Filter>(source, target, rowGrouping, choices, Filter()));
}
#endif
 
template<class T>
void reduceVectorMin(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                     std::vector<uint_fast64_t>* choices = nullptr) {
    reduceVector<T, storm::utility::ElementLess<T>>(source, target, rowGrouping, choices);
}
 
#ifdef STORM_HAVE_INTELTBB
template<class T>
void reduceVectorMinParallel(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                             std::vector<uint_fast64_t>* choices = nullptr) {
    reduceVector<T, storm::utility::ElementLess<T>>(source, target, rowGrouping, choices);
}
#endif
 
template<class T>
void reduceVectorMax(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                     std::vector<uint_fast64_t>* choices = nullptr) {
    reduceVector<T, storm::utility::ElementGreater<T>>(source, target, rowGrouping, choices);
}
 
#ifdef STORM_HAVE_INTELTBB
template<class T>
void reduceVectorMaxParallel(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping,
                             std::vector<uint_fast64_t>* choices = nullptr) {
    reduceVector<T, storm::utility::ElementGreater<T>>(source, target, rowGrouping, choices);
}
#endif
 
template<class T>
void reduceVectorMinOrMax(storm::solver::OptimizationDirection dir, std::vector<T> const& source, std::vector<T>& target,
                          std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
    if (dir == storm::solver::OptimizationDirection::Minimize) {
        reduceVectorMin(source, target, rowGrouping, choices);
    } else {
        reduceVectorMax(source, target, rowGrouping, choices);
    }
}
 
#ifdef STORM_HAVE_INTELTBB
template<class T>
void reduceVectorMinOrMaxParallel(storm::solver::OptimizationDirection dir, std::vector<T> const& source, std::vector<T>& target,
                                  std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
    if (dir == storm::solver::OptimizationDirection::Minimize) {
        reduceVectorMinParallel(source, target, rowGrouping, choices);
    } else {
        reduceVectorMaxParallel(source, target, rowGrouping, choices);
    }
}
#endif
 
template<class T>
bool equalModuloPrecision(T const& val1, T const& val2, T const& precision, bool relativeError = true) {
    if (relativeError) {
        if (storm::utility::isZero<T>(val1)) {
            return storm::utility::isZero(val2);
        }
        T relDiff = (val1 - val2) / val1;
        if (storm::utility::abs(relDiff) > precision) {
            return false;
        }
    } else {
        T diff = val1 - val2;
        if (storm::utility::abs(diff) > precision) {
            return false;
        }
    }
    return true;
}
 
// Specializiation for double as the relative check for doubles very close to zero is not meaningful.
template<>
inline bool equalModuloPrecision(double const& val1, double const& val2, double const& precision, bool relativeError) {
    if (relativeError) {
        if (storm::utility::isAlmostZero(val2)) {
            return storm::utility::isAlmostZero(val1);
        }
        double relDiff = (val1 - val2) / val1;
        if (storm::utility::abs(relDiff) > precision) {
            return false;
        }
    } else {
        double diff = val1 - val2;
        if (storm::utility::abs(diff) > precision) {
            return false;
        }
    }
    return true;
}
 
template<class T>
bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, T const& precision, bool relativeError) {
    STORM_LOG_ASSERT(vectorLeft.size() == vectorRight.size(), "Lengths of vectors does not match.");
 
    auto leftIt = vectorLeft.begin();
    auto leftIte = vectorLeft.end();
    auto rightIt = vectorRight.begin();
    for (; leftIt != leftIte; ++leftIt, ++rightIt) {
        if (!equalModuloPrecision(*leftIt, *rightIt, precision, relativeError)) {
            return false;
        }
    }
 
    return true;
}
 
template<class T>
bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, storm::storage::BitVector const& positions, T const& precision,
                          bool relativeError) {
    STORM_LOG_ASSERT(vectorLeft.size() == vectorRight.size(), "Lengths of vectors does not match.");
 
    for (auto position : positions) {
        if (!equalModuloPrecision(vectorLeft[position], vectorRight[position], precision, relativeError)) {
            return false;
        }
    }
 
    return true;
}
 
template<class T>
bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, std::vector<uint_fast64_t> const& positions, T const& precision,
                          bool relativeError) {
    STORM_LOG_ASSERT(vectorLeft.size() == vectorRight.size(), "Lengths of vectors does not match.");
 
    for (uint_fast64_t position : positions) {
        if (!equalModuloPrecision(vectorLeft[position], vectorRight[position], precision, relativeError)) {
            return false;
        }
    }
 
    return true;
}
 
template<class T>
T maximumElementAbs(std::vector<T> const& vector) {
    T res = storm::utility::zero<T>();
    for (auto const& element : vector) {
        res = std::max(res, storm::utility::abs(element));
    }
    return res;
}
 
template<class T>
T maximumElementDiff(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight) {
    T maxDiff = storm::utility::zero<T>();
    auto leftIt = vectorLeft.begin();
    auto leftIte = vectorLeft.end();
    auto rightIt = vectorRight.begin();
    for (; leftIt != leftIte; ++leftIt, ++rightIt) {
        T diff = *leftIt - *rightIt;
        T possDiff = storm::utility::abs(diff);
        maxDiff = maxDiff < possDiff ? possDiff : maxDiff;
    }
    return maxDiff;
}
 
template<class T>
T computeSquaredNorm2Difference(std::vector<T> const& b1, std::vector<T> const& b2) {
    STORM_LOG_ASSERT(b1.size() == b2.size(), "Vector sizes mismatch.");
 
    T result = storm::utility::zero<T>();
 
    auto b1It = b1.begin();
    auto b1Ite = b1.end();
    auto b2It = b2.begin();
 
    for (; b1It != b1Ite; ++b1It, ++b2It) {
        result += storm::utility::pow<T>(*b1It - *b2It, 2);
    }
 
    return result;
}
 
template<typename ValueType>
void clip(std::vector<ValueType>& x, boost::optional<ValueType> const& lowerBound, boost::optional<ValueType> const& upperBound) {
    for (auto& entry : x) {
        if (lowerBound && entry < lowerBound.get()) {
            entry = lowerBound.get();
        } else if (upperBound && entry > upperBound.get()) {
            entry = upperBound.get();
        }
    }
}
 
template<typename ValueType>
void clip(std::vector<ValueType>& x, ValueType const& bound, bool boundFromBelow) {
    for (auto& entry : x) {
        if (boundFromBelow && entry < bound) {
            entry = bound;
        } else if (!boundFromBelow && entry > bound) {
            entry = bound;
        }
    }
}
 
template<typename ValueType>
void clip(std::vector<ValueType>& x, std::vector<ValueType> const& bounds, bool boundFromBelow) {
    auto boundsIt = bounds.begin();
    for (auto& entry : x) {
        if (boundFromBelow && entry < *boundsIt) {
            entry = *boundsIt;
        } else if (!boundFromBelow && entry > *boundsIt) {
            entry = *boundsIt;
        }
        ++boundsIt;
    }
}
 
template<class T>
std::vector<T> getConstrainedOffsetVector(std::vector<T> const& offsetVector, storm::storage::BitVector const& constraint) {
    // Reserve the known amount of slots for the resulting vector.
    std::vector<uint_fast64_t> subVector(constraint.getNumberOfSetBits() + 1);
    uint_fast64_t currentRowCount = 0;
    uint_fast64_t currentIndexCount = 1;
 
    // Set the first element as this will clearly begin at offset 0.
    subVector[0] = 0;
 
    // Loop over all states that need to be kept and copy the relative indices of the nondeterministic choices over
    // to the resulting vector.
    for (auto index : constraint) {
        subVector[currentIndexCount] = currentRowCount + offsetVector[index + 1] - offsetVector[index];
        currentRowCount += offsetVector[index + 1] - offsetVector[index];
        ++currentIndexCount;
    }
 
    // Put a sentinel element at the end.
    subVector[constraint.getNumberOfSetBits()] = currentRowCount;
 
    return subVector;
}
 
template<typename TargetType, typename SourceType>
std::vector<TargetType> convertNumericVector(std::vector<SourceType> const& oldVector) {
    std::vector<TargetType> resultVector;
    resultVector.reserve(oldVector.size());
    for (auto const& oldValue : oldVector) {
        resultVector.push_back(storm::utility::convertNumber<TargetType>(oldValue));
    }
    return resultVector;
}
 
template<typename TargetType, typename SourceType>
void convertNumericVector(std::vector<SourceType> const& inputVector, std::vector<TargetType>& targetVector) {
    assert(inputVector.size() == targetVector.size());
    applyPointwise(inputVector, targetVector, [](SourceType const& v) { return storm::utility::convertNumber<TargetType>(v); });
}
 
template<typename NewValueType, typename ValueType>
std::vector<NewValueType> toValueType(std::vector<ValueType> const& oldVector) {
    std::vector<NewValueType> resultVector;
    resultVector.reserve(oldVector.size());
    for (auto const& oldValue : oldVector) {
        resultVector.push_back(static_cast<NewValueType>(oldValue));
    }
    return resultVector;
}
 
template<typename ValueType, typename TargetValueType>
typename std::enable_if<std::is_same<ValueType, storm::RationalNumber>::value, std::pair<std::vector<TargetValueType>, ValueType>>::type toIntegralVector(
    std::vector<ValueType> const& vec) {
    // Collect the numbers occurring in the input vector
    std::set<ValueType> occurringNonZeroNumbers;
    for (auto const& v : vec) {
        if (!storm::utility::isZero(v)) {
            occurringNonZeroNumbers.insert(v);
        }
    }
 
    // Compute the scaling factor
    ValueType factor;
    if (occurringNonZeroNumbers.empty()) {
        factor = storm::utility::one<ValueType>();
    } else if (occurringNonZeroNumbers.size() == 1) {
        factor = *occurringNonZeroNumbers.begin();
    } else {
        // Obtain the least common multiple of the denominators of the occurring numbers.
        // We can then multiply the numbers with the lcm to obtain integers.
        auto numberIt = occurringNonZeroNumbers.begin();
        ValueType lcm = storm::utility::asFraction(*numberIt).second;
        for (++numberIt; numberIt != occurringNonZeroNumbers.end(); ++numberIt) {
            lcm = carl::lcm(lcm, storm::utility::asFraction(*numberIt).second);
        }
        // Multiply all values with the lcm. To reduce the range of considered integers, we also obtain the gcd of the results.
        numberIt = occurringNonZeroNumbers.begin();
        ValueType gcd = *numberIt * lcm;
        for (++numberIt; numberIt != occurringNonZeroNumbers.end(); ++numberIt) {
            gcd = carl::gcd(gcd, static_cast<ValueType>(*numberIt * lcm));
        }
 
        factor = gcd / lcm;
    }
 
    // Build the result
    std::vector<TargetValueType> result;
    result.reserve(vec.size());
    for (auto const& v : vec) {
        ValueType vScaled = v / factor;
        STORM_LOG_ASSERT(storm::utility::isInteger(vScaled), "Resulting number '(" << v << ")/(" << factor << ") = " << vScaled << "' is not integral.");
        result.push_back(storm::utility::convertNumber<TargetValueType, ValueType>(vScaled));
    }
    return std::make_pair(std::move(result), std::move(factor));
}
 
template<typename ValueType, typename TargetValueType>
typename std::enable_if<!std::is_same<ValueType, storm::RationalNumber>::value, std::pair<std::vector<TargetValueType>, ValueType>>::type toIntegralVector(
    std::vector<ValueType> const& vec) {
    // TODO: avoid converting back and forth
    auto rationalNumberVec = convertNumericVector<storm::RationalNumber>(vec);
    auto rationalNumberResult = toIntegralVector<storm::RationalNumber, TargetValueType>(rationalNumberVec);
 
    return std::make_pair(std::move(rationalNumberResult.first), storm::utility::convertNumber<ValueType>(rationalNumberResult.second));
}
 
template<typename Type>
std::vector<Type> filterVector(std::vector<Type> const& in, storm::storage::BitVector const& filter) {
    std::vector<Type> result;
    result.reserve(filter.getNumberOfSetBits());
    for (auto index : filter) {
        result.push_back(in[index]);
    }
    STORM_LOG_ASSERT(result.size() == filter.getNumberOfSetBits(), "Result does not match.");
    return result;
}
 
template<typename Type>
void filterVectorInPlace(std::vector<Type>& v, storm::storage::BitVector const& filter) {
    STORM_LOG_ASSERT(v.size() == filter.size(), "The filter size does not match the size of the input vector");
    uint_fast64_t size = v.size();
    // we can start our work at the first index where the filter has value zero
    uint_fast64_t firstUnsetIndex = filter.getNextUnsetIndex(0);
    if (firstUnsetIndex < size) {
        auto vIt = v.begin() + firstUnsetIndex;
        for (uint_fast64_t index = filter.getNextSetIndex(firstUnsetIndex + 1); index != size; index = filter.getNextSetIndex(index + 1)) {
            *vIt = std::move(v[index]);
            ++vIt;
        }
        v.resize(vIt - v.begin());
        v.shrink_to_fit();
    }
    STORM_LOG_ASSERT(v.size() == filter.getNumberOfSetBits(), "Result does not match.");
}
 
template<typename T>
bool hasNegativeEntry(std::vector<T> const& v) {
    return std::any_of(v.begin(), v.end(), [](T value) { return value < storm::utility::zero<T>(); });
}
 
template<typename T>
bool hasPositiveEntry(std::vector<T> const& v) {
    return std::any_of(v.begin(), v.end(), [](T value) { return value > storm::utility::zero<T>(); });
}
 
template<typename T>
bool hasNonZeroEntry(std::vector<T> const& v) {
    return std::any_of(v.begin(), v.end(), [](T value) { return !storm::utility::isZero(value); });
}
 
template<typename T>
bool hasZeroEntry(std::vector<T> const& v) {
    return std::any_of(v.begin(), v.end(), [](T value) { return storm::utility::isZero(value); });
}
 
template<typename T>
bool hasInfinityEntry(std::vector<T> const& v) {
    return std::any_of(v.begin(), v.end(), [](T value) { return storm::utility::isInfinity(value); });
}
 
template<typename T>
std::vector<T> applyInversePermutation(std::vector<uint64_t> const& inversePermutation, std::vector<T> const& source) {
    std::vector<T> result;
    result.reserve(source.size());
    for (uint64_t sourceIndex : inversePermutation) {
        result.push_back(source[sourceIndex]);
    }
    return result;
}
 
template<typename T>
std::vector<T> applyInversePermutationToGroupedVector(std::vector<uint64_t> const& inversePermutation, std::vector<T> const& source,
                                                      std::vector<uint64_t> const& rowGroupIndices) {
    STORM_LOG_ASSERT(inversePermutation.size() == rowGroupIndices.size() - 1, "Inverse permutation and row group indices do not match.");
    std::vector<T> result;
    result.reserve(source.size());
    for (auto sourceGroupIndex : inversePermutation) {
        for (uint64_t sourceIndex = rowGroupIndices[sourceGroupIndex]; sourceIndex < rowGroupIndices[sourceGroupIndex + 1]; ++sourceIndex) {
            result.push_back(source[sourceIndex]);
        }
    }
    STORM_LOG_ASSERT(result.size() == source.size(), "result has unexpected length.");
    return result;
}
 
template<typename ValueType>
std::string toString(std::vector<ValueType> const& vector) {
    std::stringstream stream;
    stream << "vector (" << vector.size() << ") [ ";
    if (!vector.empty()) {
        for (uint_fast64_t i = 0; i < vector.size() - 1; ++i) {
            stream << vector[i] << ", ";
        }
        stream << vector.back();
    }
    stream << " ]";
    return stream.str();
}
 
template<typename PT1, typename PT2>
std::string toString(std::vector<std::pair<PT1, PT2>> const& vector) {
    std::stringstream stream;
    stream << "vector (" << vector.size() << ") [ ";
    if (!vector.empty()) {
        for (uint_fast64_t i = 0; i < vector.size() - 1; ++i) {
            stream << "{" << vector[i].first << "," << vector[i].second << "}, ";
        }
        stream << "{" << vector.back().first << "," << vector.back().second << "}";
    }
    stream << " ]";
    return stream.str();
}
 
}  // namespace vector
}  // namespace utility
}  // namespace storm
 
#endif /* STORM_UTILITY_VECTOR_H_ */