SparseMatrix.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <cstdint>
#include <iosfwd>
#include <iterator>
#include <vector>
 
#include <boost/functional/hash.hpp>
#include <boost/optional.hpp>
#include <boost/range/irange.hpp>
 
#include "storm/solver/OptimizationDirection.h"
#include "storm/storage/BitVector.h"
#include "storm/storage/sparse/StateType.h"
 
#include "storm/adapters/IntelTbbAdapter.h"
#include "storm/utility/OsDetection.h"
#include "storm/utility/constants.h"
 
// Forward declaration for adapter classes.
namespace storm {
namespace adapters {
template<typename ValueType>
class GmmxxAdapter;
class EigenAdapter;
class StormAdapter;
}  // namespace adapters
}  // namespace storm
 
namespace storm {
namespace storage {
 
// Forward declare matrix class.
template<typename T>
class SparseMatrix;
 
typedef storm::storage::sparse::state_type SparseMatrixIndexType;
 
template<typename IndexType, typename ValueType>
class MatrixEntry {
   public:
    typedef IndexType index_type;
    typedef ValueType value_type;
 
    MatrixEntry(index_type column, value_type value);
 
    MatrixEntry(std::pair<index_type, value_type>&& pair);
 
    MatrixEntry() = default;
    MatrixEntry(MatrixEntry const& other) = default;
    MatrixEntry& operator=(MatrixEntry const& other) = default;
#ifndef WINDOWS
    MatrixEntry(MatrixEntry&& other) = default;
    MatrixEntry& operator=(MatrixEntry&& other) = default;
#endif
 
    index_type const& getColumn() const;
 
    void setColumn(index_type const& column);
 
    value_type const& getValue() const;
 
    void setValue(value_type const& value);
 
    std::pair<index_type, value_type> const& getColumnValuePair() const;
 
    MatrixEntry operator*(value_type factor) const;
 
    bool operator==(MatrixEntry const& other) const;
    bool operator!=(MatrixEntry const& other) const;
 
    template<typename IndexTypePrime, typename ValueTypePrime>
    friend std::ostream& operator<<(std::ostream& out, MatrixEntry<IndexTypePrime, ValueTypePrime> const& entry);
 
   private:
    // The actual matrix entry.
    std::pair<index_type, value_type> entry;
};
 
template<typename IndexType, typename ValueType>
std::size_t hash_value(MatrixEntry<IndexType, ValueType> const& matrixEntry) {
    std::size_t seed = 0;
    boost::hash_combine(seed, matrixEntry.getColumn());
    boost::hash_combine(seed, matrixEntry.getValue());
    return seed;
}
 
template<typename ValueType>
class SparseMatrixBuilder {
   public:
    typedef SparseMatrixIndexType index_type;
    typedef ValueType value_type;
 
    SparseMatrixBuilder(index_type rows = 0, index_type columns = 0, index_type entries = 0, bool forceDimensions = true, bool hasCustomRowGrouping = false,
                        index_type rowGroups = 0);
 
    SparseMatrixBuilder(SparseMatrix<ValueType>&& matrix);
 
    void addNextValue(index_type row, index_type column, value_type const& value);
 
    void newRowGroup(index_type startingRow);
 
    /*
     * Finalizes the sparse matrix to indicate that initialization process has been completed and the matrix
     * may now be used. This must be called after all entries have been added to the matrix via addNextValue.
     *
     * @param overriddenRowCount If this is set to a value that is greater than the current number of rows,
     * this will cause the finalize method to add empty rows at the end of the matrix until the given row count
     * has been matched. Note that this will *not* override the row count that has been given upon construction
     * (if any), but will only take effect if the matrix has been created without the number of rows given.
     * @param overriddenColumnCount If this is set to a value that is greater than the current number of columns,
     * this will cause the finalize method to set the number of columns to the given value. Note that this will
     * *not* override the column count that has been given upon construction (if any), but will only take effect
     * if the matrix has been created without the number of columns given. By construction, the matrix will have
     * no entries in the columns that have been added this way.
     * @param overriddenRowGroupCount If this is set to a value that is greater than the current number of row
     * groups, this will cause the method to set the number of row groups to the given value. Note that this will
     * *not* override the row group count that has been given upon construction (if any), but will only take
     * effect if the matrix has been created without the number of row groups given. By construction, the row
     * groups added this way will be empty.
     */
    SparseMatrix<value_type> build(index_type overriddenRowCount = 0, index_type overriddenColumnCount = 0, index_type overriddenRowGroupCount = 0);
 
    index_type getLastRow() const;
 
    index_type getCurrentRowGroupCount() const;
 
    index_type getLastColumn() const;
 
    void replaceColumns(std::vector<index_type> const& replacements, index_type offset);
 
    void addDiagonalEntry(index_type row, ValueType const& value);
 
   private:
    // A flag indicating whether a row count was set upon construction.
    bool initialRowCountSet;
 
    // The row count that was initially set (if any).
    index_type initialRowCount;
 
    // A flag indicating whether a column count was set upon construction.
    bool initialColumnCountSet;
 
    // The column count that was initially set (if any).
    index_type initialColumnCount;
 
    // A flag indicating whether an entry count was set upon construction.
    bool initialEntryCountSet;
 
    // The number of entries in the matrix.
    index_type initialEntryCount;
 
    // A flag indicating whether the initially given dimensions are to be enforced on the resulting matrix.
    bool forceInitialDimensions;
 
    // A flag indicating whether the builder is to construct a custom row grouping for the matrix.
    bool hasCustomRowGrouping;
 
    // A flag indicating whether the number of row groups was set upon construction.
    bool initialRowGroupCountSet;
 
    // The number of row groups in the matrix.
    index_type initialRowGroupCount;
 
    // The vector that stores the row-group indices (if they are non-trivial).
    boost::optional<std::vector<index_type>> rowGroupIndices;
 
    // The storage for the columns and values of all entries in the matrix.
    std::vector<MatrixEntry<index_type, value_type>> columnsAndValues;
 
    // A vector containing the indices at which each given row begins. This index is to be interpreted as an
    // index in the valueStorage and the columnIndications vectors. Put differently, the values of the entries
    // in row i are valueStorage[rowIndications[i]] to valueStorage[rowIndications[i + 1]] where the last
    // entry is not included anymore.
    std::vector<index_type> rowIndications;
 
    // Stores the current number of entries in the matrix. This is used for inserting an entry into a matrix
    // with preallocated storage.
    index_type currentEntryCount;
 
    // Stores the row of the last entry in the matrix. This is used for correctly inserting an entry into a
    // matrix.
    index_type lastRow;
 
    // Stores the column of the currently last entry in the matrix. This is used for correctly inserting an
    // entry into a matrix.
    index_type lastColumn;
 
    // Stores the highest column at which an entry was inserted into the matrix.
    index_type highestColumn;
 
    // Stores the currently active row group. This is used for correctly constructing the row grouping of the
    // matrix.
    index_type currentRowGroupCount;
 
    boost::optional<ValueType> pendingDiagonalEntry;
};
 
template<typename ValueType>
class SparseMatrix {
   public:
    // Declare adapter classes as friends to use internal data.
    friend class storm::adapters::GmmxxAdapter<ValueType>;
    friend class storm::adapters::EigenAdapter;
    friend class storm::adapters::StormAdapter;
    friend class SparseMatrixBuilder<ValueType>;
 
    typedef SparseMatrixIndexType index_type;
    typedef ValueType value_type;
    typedef typename std::vector<MatrixEntry<index_type, value_type>>::iterator iterator;
    typedef typename std::vector<MatrixEntry<index_type, value_type>>::const_iterator const_iterator;
 
    class rows {
       public:
        rows(iterator begin, index_type entryCount);
 
        iterator begin();
 
        iterator end();
 
        index_type getNumberOfEntries() const;
 
       private:
        // The pointer to the columnd and value of the first entry.
        iterator beginIterator;
 
        // The number of non-zero entries in the rows.
        index_type entryCount;
    };
 
    class const_rows {
       public:
        const_rows(const_iterator begin, index_type entryCount);
 
        const_iterator begin() const;
 
        const_iterator end() const;
 
        index_type getNumberOfEntries() const;
 
       private:
        // The pointer to the column and value of the first entry.
        const_iterator beginIterator;
 
        // The number of non-zero entries in the rows.
        index_type entryCount;
    };
 
    enum MatrixStatus { UNINITIALIZED, INITIALIZED };
 
    SparseMatrix();
 
    SparseMatrix(SparseMatrix<value_type> const& other);
 
    SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements);
 
    SparseMatrix(SparseMatrix<value_type>&& other);
 
    SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications,
                 std::vector<MatrixEntry<index_type, value_type>> const& columnsAndValues, boost::optional<std::vector<index_type>> const& rowGroupIndices);
 
    SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, value_type>>&& columnsAndValues,
                 boost::optional<std::vector<index_type>>&& rowGroupIndices);
 
    SparseMatrix<value_type>& operator=(SparseMatrix<value_type> const& other);
 
    SparseMatrix<value_type>& operator=(SparseMatrix<value_type>&& other);
 
    bool operator==(SparseMatrix<value_type> const& other) const;
 
    index_type getRowCount() const;
 
    index_type getColumnCount() const;
 
    index_type getEntryCount() const;
 
    index_type getRowGroupEntryCount(index_type const group) const;
 
    index_type getNonzeroEntryCount() const;
 
    void updateNonzeroEntryCount() const;
 
    void updateDimensions() const;
 
    void updateNonzeroEntryCount(std::make_signed<index_type>::type difference);
 
    index_type getRowGroupCount() const;
 
    index_type getRowGroupSize(index_type group) const;
 
    index_type getSizeOfLargestRowGroup() const;
 
    index_type getNumRowsInRowGroups(storm::storage::BitVector const& groupConstraint) const;
 
    std::vector<index_type> const& getRowGroupIndices() const;
 
    boost::integer_range<index_type> getRowGroupIndices(index_type group) const;
 
    std::vector<index_type> swapRowGroupIndices(std::vector<index_type>&& newRowGrouping);
 
    void setRowGroupIndices(std::vector<index_type> const& newRowGroupIndices);
 
    bool hasTrivialRowGrouping() const;
 
    void makeRowGroupingTrivial();
 
    storm::storage::BitVector getRowFilter(storm::storage::BitVector const& groupConstraint) const;
 
    storm::storage::BitVector getRowFilter(storm::storage::BitVector const& groupConstraint, storm::storage::BitVector const& columnConstraints) const;
 
    storm::storage::BitVector getRowGroupFilter(storm::storage::BitVector const& rowConstraint, bool setIfForAllRowsInGroup) const;
 
    void makeRowsAbsorbing(storm::storage::BitVector const& rows, bool dropZeroEntries = false);
 
    void makeRowGroupsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, bool dropZeroEntries = false);
 
    void makeRowDirac(index_type row, index_type column, bool dropZeroEntries = false);
 
    std::vector<ValueType> getRowSumVector() const;
 
    value_type getConstrainedRowSum(index_type row, storm::storage::BitVector const& columns) const;
 
    std::vector<value_type> getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const;
 
    std::vector<value_type> getConstrainedRowGroupSumVector(storm::storage::BitVector const& rowGroupConstraint,
                                                            storm::storage::BitVector const& columnConstraint) const;
 
    SparseMatrix getSubmatrix(bool useGroups, storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint,
                              bool insertDiagonalEntries = false, storm::storage::BitVector const& makeZeroColumns = storm::storage::BitVector()) const;
 
    SparseMatrix restrictRows(storm::storage::BitVector const& rowsToKeep, bool allowEmptyRowGroups = false) const;
 
    SparseMatrix permuteRows(std::vector<index_type> const& inversePermutation) const;
 
    SparseMatrix permuteRowGroupsAndColumns(std::vector<index_type> const& inverseRowGroupPermutation, std::vector<index_type> const& columnPermutation) const;
 
    SparseMatrix filterEntries(storm::storage::BitVector const& rowFilter) const;
 
    void dropZeroEntries();
 
    bool compareRows(index_type i1, index_type i2) const;
 
    BitVector duplicateRowsInRowgroups() const;
 
    void swapRows(index_type const& row1, index_type const& row2);
 
    SparseMatrix selectRowsFromRowGroups(std::vector<index_type> const& rowGroupToRowIndexMapping, bool insertDiagonalEntries = true) const;
 
    SparseMatrix selectRowsFromRowIndexSequence(std::vector<index_type> const& rowIndexSequence, bool insertDiagonalEntries = true) const;
 
    storm::storage::SparseMatrix<value_type> transpose(bool joinGroups = false, bool keepZeros = false) const;
 
    SparseMatrix<ValueType> transposeSelectedRowsFromRowGroups(std::vector<uint64_t> const& rowGroupChoices, bool keepZeros = false) const;
 
    void convertToEquationSystem();
 
    void invertDiagonal();
 
    void negateAllNonDiagonalEntries();
 
    void deleteDiagonalEntries(bool dropZeroEntries = false);
 
    std::pair<storm::storage::SparseMatrix<value_type>, std::vector<value_type>> getJacobiDecomposition() const;
 
    template<typename OtherValueType, typename ResultValueType = OtherValueType>
    ResultValueType getPointwiseProductRowSum(storm::storage::SparseMatrix<OtherValueType> const& otherMatrix, index_type const& row) const;
 
    template<typename OtherValueType, typename ResultValueType = OtherValueType>
    std::vector<ResultValueType> getPointwiseProductRowSumVector(storm::storage::SparseMatrix<OtherValueType> const& otherMatrix) const;
 
    void multiplyWithVector(std::vector<value_type> const& vector, std::vector<value_type>& result, std::vector<value_type> const* summand = nullptr) const;
 
    void multiplyWithVectorForward(std::vector<value_type> const& vector, std::vector<value_type>& result,
                                   std::vector<value_type> const* summand = nullptr) const;
    void multiplyWithVectorBackward(std::vector<value_type> const& vector, std::vector<value_type>& result,
                                    std::vector<value_type> const* summand = nullptr) const;
#ifdef STORM_HAVE_INTELTBB
    void multiplyWithVectorParallel(std::vector<value_type> const& vector, std::vector<value_type>& result,
                                    std::vector<value_type> const* summand = nullptr) const;
#endif
 
    void multiplyAndReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector,
                           std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const;
 
    void multiplyAndReduceForward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices,
                                  std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result,
                                  std::vector<uint64_t>* choices) const;
    template<typename Compare>
    void multiplyAndReduceForward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand,
                                  std::vector<ValueType>& result, std::vector<uint64_t>* choices) const;
 
    void multiplyAndReduceBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices,
                                   std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result,
                                   std::vector<uint64_t>* choices) const;
    template<typename Compare>
    void multiplyAndReduceBackward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b,
                                   std::vector<ValueType>& result, std::vector<uint64_t>* choices) const;
#ifdef STORM_HAVE_INTELTBB
    void multiplyAndReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices,
                                   std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result,
                                   std::vector<uint64_t>* choices) const;
    template<typename Compare>
    void multiplyAndReduceParallel(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b,
                                   std::vector<ValueType>& result, std::vector<uint64_t>* choices) const;
#endif
 
    value_type multiplyRowWithVector(index_type row, std::vector<value_type> const& vector) const;
 
    void multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
 
    void scaleRowsInPlace(std::vector<value_type> const& factors);
 
    void divideRowsInPlace(std::vector<value_type> const& divisors);
 
    void performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
 
    void performWalkerChaeStep(std::vector<ValueType> const& x, std::vector<ValueType> const& columnSums, std::vector<ValueType> const& b,
                               std::vector<ValueType> const& ax, std::vector<ValueType>& result) const;
 
    value_type getRowSum(index_type row) const;
 
    index_type getNonconstantEntryCount() const;
 
    index_type getNonconstantRowGroupCount() const;
 
    bool isProbabilistic() const;
 
    bool hasOnlyPositiveEntries() const;
 
    template<typename OtherValueType>
    bool isSubmatrixOf(SparseMatrix<OtherValueType> const& matrix) const;
 
    // Returns true if the matrix is the identity matrix
    bool isIdentityMatrix() const;
 
    template<typename TPrime>
    friend std::ostream& operator<<(std::ostream& out, SparseMatrix<TPrime> const& matrix);
 
    std::string getDimensionsAsString() const;
 
    void printAsMatlabMatrix(std::ostream& out) const;
 
    std::size_t hash() const;
 
    const_rows getRows(index_type startRow, index_type endRow) const;
 
    rows getRows(index_type startRow, index_type endRow);
 
    const_rows getRow(index_type row) const;
 
    rows getRow(index_type row);
 
    const_rows getRow(index_type rowGroup, index_type offset) const;
 
    rows getRow(index_type rowGroup, index_type offset);
 
    const_rows getRowGroup(index_type rowGroup) const;
 
    rows getRowGroup(index_type rowGroup);
 
    const_iterator begin(index_type row) const;
 
    iterator begin(index_type row);
 
    const_iterator begin() const;
 
    iterator begin();
 
    const_iterator end(index_type row) const;
 
    iterator end(index_type row);
 
    const_iterator end() const;
 
    iterator end();
 
    template<typename NewValueType>
    SparseMatrix<NewValueType> toValueType() const {
        std::vector<MatrixEntry<SparseMatrix::index_type, NewValueType>> newColumnsAndValues;
        std::vector<SparseMatrix::index_type> newRowIndications(rowIndications);
        boost::optional<std::vector<SparseMatrix::index_type>> newRowGroupIndices(rowGroupIndices);
 
        newColumnsAndValues.resize(columnsAndValues.size());
        std::transform(
            columnsAndValues.begin(), columnsAndValues.end(), newColumnsAndValues.begin(), [](MatrixEntry<SparseMatrix::index_type, ValueType> const& a) {
                return MatrixEntry<SparseMatrix::index_type, NewValueType>(a.getColumn(), storm::utility::convertNumber<NewValueType, ValueType>(a.getValue()));
            });
 
        return SparseMatrix<NewValueType>(columnCount, std::move(newRowIndications), std::move(newColumnsAndValues), std::move(newRowGroupIndices));
    }
 
   private:
    SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint,
                              std::vector<index_type> const& rowGroupIndices, bool insertDiagonalEntries = false,
                              storm::storage::BitVector const& makeZeroColumns = storm::storage::BitVector()) const;
 
    // The number of rows of the matrix.
    index_type rowCount;
 
    // The number of columns of the matrix.
    mutable index_type columnCount;
 
    // The number of entries in the matrix.
    index_type entryCount;
 
    // The number of nonzero entries in the matrix.
    mutable index_type nonzeroEntryCount;
 
    // The storage for the columns and values of all entries in the matrix.
    std::vector<MatrixEntry<index_type, value_type>> columnsAndValues;
 
    // A vector containing the indices at which each given row begins. This index is to be interpreted as an
    // index in the valueStorage and the columnIndications vectors. Put differently, the values of the entries
    // in row i are valueStorage[rowIndications[i]] to valueStorage[rowIndications[i + 1]] where the last
    // entry is not included anymore.
    std::vector<index_type> rowIndications;
 
    // A flag indicating whether the matrix has a trivial row grouping. Note that this may be true and yet
    // there may be row group indices, because they were requested from the outside.
    bool trivialRowGrouping;
 
    // A vector indicating the row groups of the matrix. This needs to be mutible in case we create it on-the-fly.
    mutable boost::optional<std::vector<index_type>> rowGroupIndices;
};
 
}  // namespace storage
}  // namespace storm