sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Made changes to all files that use the Sparse Matrix, as the diagonal entries are now included in the main storage. This refs #34

main
PBerger 13 years ago
parent
11b16f5dde
commit
b89db58fbb
8 changed files with 305 additions and 359 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 45
      src/adapters/GmmxxAdapter.h
  2. 2
      src/formula/BoundedEventually.h
  3. 8
      src/models/GraphTransitions.h
  4. 10
      src/parser/DeterministicSparseTransitionParser.cpp
  5. 462
      src/storage/SquareSparseMatrix.h
  6. 113
      src/utility/IoUtility.cpp
  7. 4
      test/parser/ReadTraFileTest.cpp
  8. 20
      test/storage/SquareSparseMatrixTest.cpp

45
src/adapters/GmmxxAdapter.h
View File

@ -27,7 +27,7 @@ public:
*/
template<class T>
static gmm::csr_matrix<T>* toGmmxxSparseMatrix(storm::storage::SquareSparseMatrix<T> const& matrix) {
uint_fast64_t realNonZeros = matrix.getNonZeroEntryCount() + matrix.getDiagonalNonZeroEntryCount();
uint_fast64_t realNonZeros = matrix.getNonZeroEntryCount();
LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
// Prepare the resulting matrix.
@ -44,48 +44,23 @@ public:
T* tmpValueArray = new T[realNonZeros];
T zero(0);
uint_fast64_t currentPosition = 0;
uint_fast64_t insertedDiagonalElements = 0;
for (uint_fast64_t i = 0; i < matrix.rowCount; ++i) {
// Compute correct start index of row.
result->jc[i] = matrix.rowIndications[i] + insertedDiagonalElements;
// If the current row has no non-zero which is not on the diagonal, we have to check the
// diagonal element explicitly.
if (matrix.rowIndications[i + 1] - matrix.rowIndications[i] == 0) {
if (matrix.diagonalStorage[i] != zero) {
tmpColumnIndicationsArray[currentPosition] = i;
tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
++currentPosition; ++insertedDiagonalElements;
}
} else {
// Otherwise, we can just enumerate the non-zeros which are not on the diagonal
// and fit in the diagonal element where appropriate.
bool includedDiagonal = false;
for (uint_fast64_t j = matrix.rowIndications[i]; j < matrix.rowIndications[i + 1]; ++j) {
if (matrix.diagonalStorage[i] != zero && !includedDiagonal && matrix.columnIndications[j] > i) {
includedDiagonal = true;
tmpColumnIndicationsArray[currentPosition] = i;
tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
++currentPosition; ++insertedDiagonalElements;
}
tmpColumnIndicationsArray[currentPosition] = matrix.columnIndications[j];
tmpValueArray[currentPosition] = matrix.valueStorage[j];
++currentPosition;
}
// If the diagonal element is non-zero and was not inserted until now (i.e. all
// off-diagonal elements in the row are before the diagonal element.
if (!includedDiagonal && matrix.diagonalStorage[i] != zero) {
tmpColumnIndicationsArray[currentPosition] = i;
tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
++currentPosition; ++insertedDiagonalElements;
}
result->jc[i] = matrix.rowIndications[i];
// Otherwise, we can just enumerate the non-zeros which are not on the diagonal
// and fit in the diagonal element where appropriate.
for (uint_fast64_t j = matrix.rowIndications[i]; j < matrix.rowIndications[i + 1]; ++j) {
tmpColumnIndicationsArray[currentPosition] = matrix.columnIndications[j];
tmpValueArray[currentPosition] = matrix.valueStorage[j];
++currentPosition;
}
}
// Fill in sentinel element at the end.
result->jc[matrix.rowCount] = static_cast<unsigned int>(realNonZeros);
// Now, we can copy the temporary array to the GMMXX format.
// FIXME Now everything can just be copied. No TMP needed anymore
result->ir.resize(realNonZeros);
std::copy(tmpColumnIndicationsArray, tmpColumnIndicationsArray + realNonZeros, result->ir.begin());
delete[] tmpColumnIndicationsArray;

2
src/formula/BoundedEventually.h
View File

@ -120,7 +120,7 @@ public:
BoundedEventually<T>* result = new BoundedEventually<T>();
result->setBound(bound);
if (child != nullptr) {
result->setRight(child->clone());
result->setChild(child->clone());
}
return result;
}

8
src/models/GraphTransitions.h
View File

@ -93,12 +93,12 @@ private:
// First, we copy the index list from the sparse matrix as this will
// stay the same.
std::copy(transitionMatrix->getRowIndicationsPointer(), transitionMatrix->getRowIndicationsPointer() + numberOfStates + 1, this->stateIndications);
std::copy(transitionMatrix->getRowIndicationsPointer().begin(), transitionMatrix->getRowIndicationsPointer().end(), this->stateIndications);
// Now we can iterate over all rows of the transition matrix and record
// the target state.
for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) {
for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
this->stateIndications[currentNonZeroElement++] = *rowIt;
}
}
@ -117,7 +117,7 @@ private:
// NOTE: We disregard the diagonal here, as we only consider "true"
// predecessors.
for (uint_fast64_t i = 0; i < numberOfStates; i++) {
for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
this->stateIndications[*rowIt + 1]++;
}
}
@ -140,7 +140,7 @@ private:
// Now we are ready to actually fill in the list of predecessors for
// every state. Again, we start by considering all but the last row.
for (uint_fast64_t i = 0; i < numberOfStates; i++) {
for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
this->successorList[nextIndicesList[*rowIt]++] = i;
}
}

10
src/parser/DeterministicSparseTransitionParser.cpp
View File

@ -32,13 +32,10 @@ namespace parser{
* non-zero cells and maximum node id.
*
* This method does the first pass through the .tra file and computes
* the number of non-zero elements that are not diagonal elements,
* which correspondents to the number of transitions that are not
* self-loops.
* (Diagonal elements are treated in a special way).
* the number of non-zero elements.
* It also calculates the maximum node id and stores it in maxnode.
*
* @return The number of non-zero elements that are not on the diagonal
* @return The number of non-zero elements
* @param buf Data to scan. Is expected to be some char array.
* @param maxnode Is set to highest id of all nodes.
*/
@ -89,7 +86,8 @@ uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, uint_fas
LOG4CPLUS_ERROR(logger, "Expected a positive probability but got \"" << std::string(buf, 0, 16) << "\".");
return 0;
}
if (row == col) non_zero--;
// not needed anymore
//if (row == col) non_zero--;
buf = trimWhitespaces(tmp);
}

462
src/storage/SquareSparseMatrix.h
View File

@ -5,6 +5,7 @@
#include <new>
#include <algorithm>
#include <iostream>
#include <iterator>
#include "boost/integer/integer_mask.hpp"
#include "src/exceptions/InvalidStateException.h"
@ -31,8 +32,7 @@ namespace storm {
namespace storage {
/*!
* A sparse matrix class with a constant number of non-zero entries on the non-diagonal fields
* and a separate dense storage for the diagonal elements.
* A sparse matrix class with a constant number of non-zero entries.
* NOTE: Addressing *is* zero-based, so the valid range for getValue and addNextValue is 0..(rows - 1)
* where rows is the first argument to the constructor.
*/
@ -73,9 +73,25 @@ public:
* Constructs a sparse matrix object with the given number of rows.
* @param rows The number of rows of the matrix
*/
SquareSparseMatrix(uint_fast64_t rows)
: rowCount(rows), nonZeroEntryCount(0), valueStorage(nullptr),
diagonalStorage(nullptr),columnIndications(nullptr), rowIndications(nullptr),
SquareSparseMatrix(uint_fast64_t rows, uint_fast64_t cols)
: rowCount(rows), colCount(cols), nonZeroEntryCount(0),
internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
/* Sadly, Delegate Constructors are not yet available with MSVC2012
//! Constructor
/*!
* Constructs a square sparse matrix object with the given number rows
* @param size The number of rows and cols in the matrix
*/ /*
SquareSparseMatrix(uint_fast64_t size) : SquareSparseMatrix(size, size) { }
*/
//! Constructor
/*!
* Constructs a square sparse matrix object with the given number rows
* @param size The number of rows and cols in the matrix
*/
SquareSparseMatrix(uint_fast64_t size) : rowCount(size), colCount(size), nonZeroEntryCount(0),
internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
//! Copy Constructor
@ -84,7 +100,7 @@ public:
* @param ssm A reference to the matrix to be copied.
*/
SquareSparseMatrix(const SquareSparseMatrix<T> &ssm)
: rowCount(ssm.rowCount), nonZeroEntryCount(ssm.nonZeroEntryCount),
: rowCount(ssm.rowCount), colCount(ssm.colCount), nonZeroEntryCount(ssm.nonZeroEntryCount),
internalStatus(ssm.internalStatus), currentSize(ssm.currentSize), lastRow(ssm.lastRow) {
LOG4CPLUS_WARN(logger, "Invoking copy constructor.");
// Check whether copying the matrix is safe.
@ -98,24 +114,12 @@ public:
LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
throw std::bad_alloc();
} else {
// Now that all storages have been prepared, copy over all
// elements. Start by copying the elements of type value and
// copy them seperately in order to invoke copy their copy
// constructor. This may not be necessary, but it is safer to
// do so in any case.
for (uint_fast64_t i = 0; i < nonZeroEntryCount; ++i) {
// use T() to force use of the copy constructor for complex T types
valueStorage[i] = T(ssm.valueStorage[i]);
}
for (uint_fast64_t i = 0; i < rowCount; ++i) {
// use T() to force use of the copy constructor for complex T types
diagonalStorage[i] = T(ssm.diagonalStorage[i]);
}
std::copy(ssm.valueStorage.begin(), ssm.valueStorage.end(), std::back_inserter(valueStorage));
// The elements that are not of the value type but rather the
// index type may be copied directly.
std::copy(ssm.columnIndications, ssm.columnIndications + nonZeroEntryCount, columnIndications);
std::copy(ssm.rowIndications, ssm.rowIndications + rowCount + 1, rowIndications);
std::copy(ssm.columnIndications.begin(), ssm.columnIndications.end(), std::back_inserter(columnIndications));
std::copy(ssm.rowIndications.begin(), ssm.rowIndications.end(), std::back_inserter(rowIndications));
}
}
}
@ -126,18 +130,9 @@ public:
*/
~SquareSparseMatrix() {
setState(MatrixStatus::UnInitialized);
if (valueStorage != nullptr) {
delete[] valueStorage;
}
if (columnIndications != nullptr) {
delete[] columnIndications;
}
if (rowIndications != nullptr) {
delete[] rowIndications;
}
if (diagonalStorage != nullptr) {
delete[] diagonalStorage;
}
valueStorage.resize(0);
columnIndications.resize(0);
rowIndications.resize(0);
}
/*!
@ -155,11 +150,11 @@ public:
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to initialize matrix that is not uninitialized.");
throw storm::exceptions::InvalidStateException("Trying to initialize matrix that is not uninitialized.");
} else if (rowCount == 0) {
} else if ((rowCount == 0) || (colCount == 0)) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to create initialize a matrix with 0 rows.");
throw storm::exceptions::InvalidArgumentException("Trying to create initialize a matrix with 0 rows.");
} else if (((rowCount * rowCount) - rowCount) < nonZeroEntries) {
LOG4CPLUS_ERROR(logger, "Trying to create initialize a matrix with 0 rows or 0 columns.");
throw storm::exceptions::InvalidArgumentException("Trying to create initialize a matrix with 0 rows or 0 columns.");
} else if ((rowCount * colCount) < nonZeroEntries) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to initialize a matrix with more non-zero entries than there can be.");
throw storm::exceptions::InvalidArgumentException("Trying to initialize a matrix with more non-zero entries than there can be.");
@ -196,37 +191,72 @@ public:
throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that is not in compressed form.");
}
// Compute the actual (i.e. non-diagonal) number of non-zero entries.
nonZeroEntryCount = getEigenSparseMatrixCorrectNonZeroEntryCount(eigenSparseMatrix);
if (eigenSparseMatrix.rows() > this->rowCount) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to initialize from an Eigen matrix that has more rows than the target matrix.");
throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that has more rows than the target matrix.");
}
if (eigenSparseMatrix.cols() > this->colCount) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to initialize from an Eigen matrix that has more columns than the target matrix.");
throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that has more columns than the target matrix.");
}
const _Index entryCount = eigenSparseMatrix.nonZeros();
nonZeroEntryCount = entryCount;
lastRow = 0;
// Try to prepare the internal storage and throw an error in case of
// failure.
if (!prepareInternalStorage()) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
throw std::bad_alloc();
} else {
// Get necessary pointers to the contents of the Eigen matrix.
const T* valuePtr = eigenSparseMatrix.valuePtr();
const _Index* indexPtr = eigenSparseMatrix.innerIndexPtr();
const _Index* outerPtr = eigenSparseMatrix.outerIndexPtr();
// If the given matrix is in RowMajor format, copying can simply
// be done by adding all values in order.
// Direct copying is, however, prevented because we have to
// separate the diagonal entries from others.
if (isEigenRowMajor(eigenSparseMatrix)) {
// Because of the RowMajor format outerSize evaluates to the
// number of rows.
const _Index rowCount = eigenSparseMatrix.outerSize();
for (_Index row = 0; row < rowCount; ++row) {
for (_Index col = outerPtr[row]; col < outerPtr[row + 1]; ++col) {
addNextValue(row, indexPtr[col], valuePtr[col]);
}
// Get necessary pointers to the contents of the Eigen matrix.
const T* valuePtr = eigenSparseMatrix.valuePtr();
const _Index* indexPtr = eigenSparseMatrix.innerIndexPtr();
const _Index* outerPtr = eigenSparseMatrix.outerIndexPtr();
// If the given matrix is in RowMajor format, copying can simply
// be done by adding all values in order.
// Direct copying is, however, prevented because we have to
// separate the diagonal entries from others.
if (isEigenRowMajor(eigenSparseMatrix)) {
// Because of the RowMajor format outerSize evaluates to the
// number of rows.
if (!prepareInternalStorage(false)) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
throw std::bad_alloc();
} else {
if ((eigenSparseMatrix.innerSize() > nonZeroEntryCount) || (entryCount > nonZeroEntryCount)) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Invalid internal composition of Eigen Sparse Matrix");
throw storm::exceptions::InvalidArgumentException("Invalid internal composition of Eigen Sparse Matrix");
}
std::vector<uint_fast64_t> eigenColumnTemp;
std::vector<uint_fast64_t> eigenRowTemp;
std::vector<T> eigenValueTemp;
uint_fast64_t outerSize = eigenSparseMatrix.outerSize() + 1;
for (uint_fast64_t i = 0; i < entryCount; ++i) {
eigenColumnTemp.push_back(indexPtr[i]);
eigenValueTemp.push_back(valuePtr[i]);
}
for (uint_fast64_t i = 0; i < outerSize; ++i) {
eigenRowTemp.push_back(outerPtr[i]);
}
std::copy(eigenRowTemp.begin(), eigenRowTemp.end(), std::back_inserter(this->rowIndications));
std::copy(eigenColumnTemp.begin(), eigenColumnTemp.end(), std::back_inserter(this->columnIndications));
std::copy(eigenValueTemp.begin(), eigenValueTemp.end(), std::back_inserter(this->valueStorage));
currentSize = entryCount;
lastRow = rowCount;
}
} else {
if (!prepareInternalStorage()) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
throw std::bad_alloc();
} else {
const _Index entryCount = eigenSparseMatrix.nonZeros();
// Because of the ColMajor format outerSize evaluates to the
// number of columns.
const _Index colCount = eigenSparseMatrix.outerSize();
@ -250,8 +280,7 @@ public:
// add it in case it is also in the current row.
if ((positions[currentColumn] < outerPtr[currentColumn + 1])
&& (indexPtr[positions[currentColumn]] == currentRow)) {
addNextValue(currentRow, currentColumn,
valuePtr[positions[currentColumn]]);
addNextValue(currentRow, currentColumn, valuePtr[positions[currentColumn]]);
// Remember that we found one more non-zero element.
++i;
// Mark this position as "used".
@ -268,8 +297,8 @@ public:
}
delete[] positions;
}
setState(MatrixStatus::Initialized);
}
setState(MatrixStatus::Initialized);
}
/*!
@ -283,30 +312,27 @@ public:
void addNextValue(const uint_fast64_t row, const uint_fast64_t col, const T& value) {
// Check whether the given row and column positions are valid and throw
// error otherwise.
if ((row > rowCount) || (col > rowCount)) {
if ((row > rowCount) || (col > colCount)) {
triggerErrorState();
LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ").");
throw storm::exceptions::OutOfRangeException("Trying to add a value at illegal position.");
}
if (row == col) { // Set a diagonal element.
diagonalStorage[row] = value;
} else { // Set a non-diagonal element.
// If we switched to another row, we have to adjust the missing
// entries in the row_indications array.
if (row != lastRow) {
for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
rowIndications[i] = currentSize;
}
lastRow = row;
// If we switched to another row, we have to adjust the missing
// entries in the row_indications array.
if (row != lastRow) {
for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
rowIndications[i] = currentSize;
}
lastRow = row;
}
// Finally, set the element and increase the current size.
valueStorage[currentSize] = value;
columnIndications[currentSize] = col;
// Finally, set the element and increase the current size.
valueStorage[currentSize] = value;
columnIndications[currentSize] = col;
++currentSize;
}
++currentSize;
}
/*
@ -355,18 +381,12 @@ public:
*/
inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) const {
// Check for illegal access indices.
if ((row > rowCount) || (col > rowCount)) {
if ((row > rowCount) || (col > colCount)) {
LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ").");
throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position.");
return false;
}
// Read elements on the diagonal directly.
if (row == col) {
*target = diagonalStorage[row];
return true;
}
// In case the element is not on the diagonal, we have to iterate
// over the accessed row to find the element.
uint_fast64_t rowStart = rowIndications[row];
@ -405,17 +425,12 @@ public:
*/
inline T& getValue(uint_fast64_t row, uint_fast64_t col) {
// Check for illegal access indices.
if ((row > rowCount) || (col > rowCount)) {
if ((row > rowCount) || (col > colCount)) {
LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ").");
throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position.");
}
// Read elements on the diagonal directly.
if (row == col) {
return diagonalStorage[row];
}
// In case the element is not on the diagonal, we have to iterate
// we have to iterate
// over the accessed row to find the element.
uint_fast64_t rowStart = rowIndications[row];
uint_fast64_t rowEnd = rowIndications[row + 1];
@ -445,20 +460,18 @@ public:
}
/*!
* Returns a pointer to the value storage of the matrix. This storage does
* *not* include elements on the diagonal.
* @return A pointer to the value storage of the matrix.
* Returns the number of columns of the matrix.
*/
T* getStoragePointer() const {
return valueStorage;
uint_fast64_t getColumnCount() const {
return colCount;
}
/*!
* Returns a pointer to the storage of elements on the diagonal.
* @return A pointer to the storage of elements on the diagonal.
* Returns a pointer to the value storage of the matrix.
* @return A pointer to the value storage of the matrix.
*/
T* getDiagonalStoragePointer() const {
return diagonalStorage;
std::vector<T>& getStoragePointer() {
return valueStorage;
}
/*!
@ -467,17 +480,17 @@ public:
* @return A pointer to the array that stores the start indices of non-zero
* entries in the value storage for each row.
*/
uint_fast64_t* getRowIndicationsPointer() const {
std::vector<uint_fast64_t>& getRowIndicationsPointer() {
return rowIndications;
}
/*!
* Returns a pointer to an array that stores the column of each non-zero
* element that is not on the diagonal.
* element.
* @return A pointer to an array that stores the column of each non-zero
* element that is not on the diagonal.
* element.
*/
uint_fast64_t* getColumnIndicationsPointer() const {
std::vector<uint_fast64_t>& getColumnIndicationsPointer() {
return columnIndications;
}
@ -548,10 +561,6 @@ public:
#define STORM_USE_TRIPLETCONVERT
# ifdef STORM_USE_TRIPLETCONVERT
// FIXME: Wouldn't it be more efficient to add the elements in
// order including the diagonal elements? Otherwise, Eigen has to
// perform some sorting.
// Prepare the triplet storage.
typedef Eigen::Triplet<T> IntTriplet;
std::vector<IntTriplet> tripletList;
@ -572,12 +581,6 @@ public:
}
}
// Then add the elements on the diagonal.
for (uint_fast64_t i = 0; i < rowCount; ++i) {
if (diagonalStorage[i] == 0) zeroCount++;
tripletList.push_back(IntTriplet(static_cast<int_fast32_t>(i), static_cast<int_fast32_t>(i), diagonalStorage[i]));
}
// Let Eigen create a matrix from the given list of triplets.
mat->setFromTriplets(tripletList.begin(), tripletList.end());
@ -596,10 +599,6 @@ public:
rowStart = rowIndications[row];
rowEnd = rowIndications[row + 1];
// Insert the element on the diagonal.
mat->insert(row, row) = diagonalStorage[row];
count++;
// Insert the elements that are not on the diagonal
while (rowStart < rowEnd) {
mat->insert(row, columnIndications[rowStart]) = valueStorage[rowStart];
@ -628,19 +627,6 @@ public:
return nonZeroEntryCount;
}
/*!
* Returns the number of non-zero entries on the diagonal.
* @return The number of non-zero entries on the diagonal.
*/
uint_fast64_t getDiagonalNonZeroEntryCount() const {
uint_fast64_t result = 0;
T zero(0);
for (uint_fast64_t i = 0; i < rowCount; ++i) {
if (diagonalStorage[i] != zero) ++result;
}
return result;
}
/*!
* This function makes the rows given by the bit vector absorbing.
* @param rows A bit vector indicating which rows to make absorbing.
@ -658,7 +644,7 @@ public:
/*!
* This function makes the given row absorbing. This means that all
* entries in will be set to 0 and the value 1 will be written
* to the element on the diagonal.
* to the element on the (pseudo-) diagonal.
* @param row The row to be made absorbing.
* @returns True iff the operation was successful.
*/
@ -675,13 +661,31 @@ public:
uint_fast64_t rowStart = rowIndications[row];
uint_fast64_t rowEnd = rowIndications[row + 1];
if (rowStart >= rowEnd) {
LOG4CPLUS_ERROR(logger, "The row " << row << " can not be made absorbing, no state in row, would have to recreate matrix!");
throw storm::exceptions::InvalidStateException("A row can not be made absorbing, no state in row, would have to recreate matrix!");
}
uint_fast64_t pseudoDiagonal = row % colCount;
bool foundDiagonal = false;
while (rowStart < rowEnd) {
valueStorage[rowStart] = storm::utility::constGetZero<T>();
if (!foundDiagonal && columnIndications[rowStart] >= pseudoDiagonal) {
foundDiagonal = true;
// insert/replace the diagonal entry
columnIndications[rowStart] = pseudoDiagonal;
valueStorage[rowStart] = storm::utility::constGetOne<T>();
} else {
valueStorage[rowStart] = storm::utility::constGetZero<T>();
}
++rowStart;
}
// Set the element on the diagonal to one.
diagonalStorage[row] = storm::utility::constGetOne<T>();
if (!foundDiagonal) {
--rowStart;
columnIndications[rowStart] = pseudoDiagonal;
valueStorage[rowStart] = storm::utility::constGetOne<T>();
}
return true;
}
@ -763,8 +767,6 @@ public:
// Copy over selected entries.
uint_fast64_t rowCount = 0;
for (auto rowIndex : constraint) {
result->addNextValue(rowCount, rowCount, diagonalStorage[rowIndex]);
for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
if (constraint.get(columnIndications[i])) {
result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]);
@ -794,8 +796,17 @@ public:
*/
void invertDiagonal() {
T one(1);
for (uint_fast64_t i = 0; i < rowCount; ++i) {
diagonalStorage[i] = one - diagonalStorage[i];
for (uint_fast64_t row = 0; row < rowCount; ++row) {
uint_fast64_t rowStart = rowIndications[row];
uint_fast64_t rowEnd = rowIndications[row + 1];
uint_fast64_t pseudoDiagonal = row % colCount;
while (rowStart < rowEnd) {
if (columnIndications[rowStart] == pseudoDiagonal) {
valueStorage[rowStart] = one - valueStorage[rowStart];
break;
}
++rowStart;
}
}
}
@ -803,8 +814,16 @@ public:
* Negates all non-zero elements that are not on the diagonal.
*/
void negateAllNonDiagonalElements() {
for (uint_fast64_t i = 0; i < nonZeroEntryCount; ++i) {
valueStorage[i] = - valueStorage[i];
for (uint_fast64_t row = 0; row < rowCount; ++row) {
uint_fast64_t rowStart = rowIndications[row];
uint_fast64_t rowEnd = rowIndications[row + 1];
uint_fast64_t pseudoDiagonal = row % colCount;
while (rowStart < rowEnd) {
if (columnIndications[rowStart] != pseudoDiagonal) {
valueStorage[rowStart] = - valueStorage[rowStart];
}
++rowStart;
}
}
}
@ -844,7 +863,6 @@ public:
// in case the given matrix does not have a non-zero element at this column position, or
// multiply the two entries and add the result to the corresponding position in the vector.
for (uint_fast64_t row = 0; row < rowCount && row < otherMatrix.rowCount; ++row) {
(*result)[row] += diagonalStorage[row] * otherMatrix.diagonalStorage[row];
for (uint_fast64_t element = rowIndications[row], nextOtherElement = otherMatrix.rowIndications[row]; element < rowIndications[row + 1] && nextOtherElement < otherMatrix.rowIndications[row + 1]; ++element) {
if (columnIndications[element] < otherMatrix.columnIndications[nextOtherElement]) {
continue;
@ -868,25 +886,23 @@ public:
uint_fast64_t getSizeInMemory() const {
uint_fast64_t size = sizeof(*this);
// Add value_storage size.
size += sizeof(T) * nonZeroEntryCount;
// Add diagonal_storage size.
size += sizeof(T) * (rowCount + 1);
size += sizeof(T) * valueStorage.capacity();
// Add column_indications size.
size += sizeof(uint_fast64_t) * nonZeroEntryCount;
size += sizeof(uint_fast64_t) * columnIndications.capacity();
// Add row_indications size.
size += sizeof(uint_fast64_t) * (rowCount + 1);
size += sizeof(uint_fast64_t) * rowIndications.capacity();
return size;
}
/*!
* Returns an iterator to the columns of the non-zero entries of the given
* row that are not on the diagonal.
* row.
* @param row The row whose columns the iterator will return.
* @return An iterator to the columns of the non-zero entries of the given
* row that are not on the diagonal.
* row.
*/
constIndexIterator beginConstColumnNoDiagIterator(uint_fast64_t row) const {
return this->columnIndications + this->rowIndications[row];
constIndexIterator beginConstColumnIterator(uint_fast64_t row) const {
return &(this->columnIndications[0]) + this->rowIndications[row];
}
/*!
@ -894,18 +910,18 @@ public:
* @param row The row for which the iterator should point to the past-the-end
* element.
*/
constIndexIterator endConstColumnNoDiagIterator(uint_fast64_t row) const {
return this->columnIndications + this->rowIndications[row + 1];
constIndexIterator endConstColumnIterator(uint_fast64_t row) const {
return &(this->columnIndications[0]) + this->rowIndications[row + 1];
}
/*!
* Returns an iterator over the elements of the given row. The iterator
* will include neither the diagonal element nor zero entries.
* will include no zero entries.
* @param row The row whose elements the iterator will return.
* @return An iterator over the elements of the given row.
*/
constIterator beginConstNoDiagIterator(uint_fast64_t row) const {
return this->valueStorage + this->rowIndications[row];
constIterator beginConstIterator(uint_fast64_t row) const {
return &(this->valueStorage[0]) + this->rowIndications[row];
}
/*!
* Returns an iterator pointing to the first element after the given
@ -914,32 +930,28 @@ public:
* past-the-end element.
* @return An iterator to the element after the given row.
*/
constIterator endConstNoDiagIterator(uint_fast64_t row) const {
return this->valueStorage + this->rowIndications[row + 1];
constIterator endConstIterator(uint_fast64_t row) const {
return &(this->valueStorage[0]) + this->rowIndications[row + 1];
}
/*!
* @brief Calculate sum of all entries in given row.
*
* Adds up all values in the given row (including the diagonal value)
* Adds up all values in the given row
* and returns the sum.
* @param row The row that should be added up.
* @return Sum of the row.
*/
T getRowSum(uint_fast64_t row) const {
T sum = this->diagonalStorage[row];
for (auto it = this->beginConstNoDiagIterator(row); it != this->endConstNoDiagIterator(row); it++) {
T sum = storm::utility::constGetZero<T>();
for (auto it = this->beginConstIterator(row); it != this->endConstIterator(row); it++) {
sum += *it;
}
return sum;
}
void print() const {
std::cout << "diag: --------------------------------" << std::endl;
for (uint_fast64_t i = 0; i < rowCount; ++i) {
std::cout << "(" << i << "," << i << ") = " << diagonalStorage[i] << std::endl;
}
std::cout << "non diag: ----------------------------" << std::endl;
std::cout << "entries: ----------------------------" << std::endl;
for (uint_fast64_t i = 0; i < rowCount; ++i) {
for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
std::cout << "(" << i << "," << columnIndications[j] << ") = " << valueStorage[j] << std::endl;
@ -955,31 +967,31 @@ private:
uint_fast64_t rowCount;
/*!
* The number of non-zero elements that are not on the diagonal.
* The number of columns of the matrix.
*/
uint_fast64_t nonZeroEntryCount;
uint_fast64_t colCount;
/*!
* Stores all non-zero values that are not on the diagonal.
* The number of non-zero elements.
*/
T* valueStorage;
uint_fast64_t nonZeroEntryCount;
/*!
* Stores all elements on the diagonal, even the ones that are zero.
* Stores all non-zero values.
*/
T* diagonalStorage;
std::vector<T> valueStorage;
/*!
* Stores the column for each non-zero element that is not on the diagonal.
* Stores the column for each non-zero element.
*/
uint_fast64_t* columnIndications;
std::vector<uint_fast64_t> columnIndications;
/*!
* Array containing the boundaries (indices) in the value_storage array
* Vector containing the boundaries (indices) in the value_storage array
* for each row. All elements of value_storage with indices between the
* i-th and the (i+1)-st element of this array belong to row i.
*/
uint_fast64_t* rowIndications;
std::vector<uint_fast64_t> rowIndications;
/*!
* The internal status of the matrix.
@ -1017,24 +1029,37 @@ private:
/*!
* Prepares the internal CSR storage. For this, it requires
* non_zero_entry_count and row_count to be set correctly.
* @param alsoPerformAllocation If set to true, all entries are pre-allocated. This is the default.
* @return True on success, false otherwise (allocation failed).
*/
bool prepareInternalStorage() {
// Set up the arrays for the elements that are not on the diagonal.
valueStorage = new (std::nothrow) T[nonZeroEntryCount]();
columnIndications = new (std::nothrow) uint_fast64_t[nonZeroEntryCount]();
// Set up the row_indications array and reserve one element more than
// there are rows in order to put a sentinel element at the end,
// which eases iteration process.
rowIndications = new (std::nothrow) uint_fast64_t[rowCount + 1]();
// Set up the array for the elements on the diagonal.
diagonalStorage = new (std::nothrow) T[rowCount]();
bool prepareInternalStorage(const bool alsoPerformAllocation) {
if (alsoPerformAllocation) {
// Set up the arrays for the elements that are not on the diagonal.
valueStorage.resize(nonZeroEntryCount, storm::utility::constGetZero<T>());
columnIndications.resize(nonZeroEntryCount, 0);
// Set up the row_indications vector and reserve one element more than
// there are rows in order to put a sentinel element at the end,
// which eases iteration process.
rowIndications.resize(rowCount + 1, 0);
// Return whether all the allocations could be made without error.
return ((valueStorage.capacity() >= nonZeroEntryCount) && (columnIndications.capacity() >= nonZeroEntryCount)
&& (rowIndications.capacity() >= (rowCount + 1)));
} else {
valueStorage.reserve(nonZeroEntryCount);
columnIndications.reserve(nonZeroEntryCount);
rowIndications.reserve(rowCount + 1);
return true;
}
}
// Return whether all the allocations could be made without error.
return ((valueStorage != NULL) && (columnIndications != NULL)
&& (rowIndications != NULL) && (diagonalStorage != NULL));
/*!
* Shorthand for prepareInternalStorage(true)
* @see prepareInternalStorage(const bool)
*/
bool prepareInternalStorage() {
return this->prepareInternalStorage(true);
}
/*!
@ -1060,53 +1085,6 @@ private:
return false;
}
/*!
* Helper function to determine the number of non-zero elements that are
* not on the diagonal of the given Eigen matrix.
* @param eigen_sparse_matrix The Eigen matrix to analyze.
* @return The number of non-zero elements that are not on the diagonal of
* the given Eigen matrix.
*/
template<typename _Scalar, int _Options, typename _Index>
_Index getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index>& eigen_sparse_matrix) const {
const _Index* indexPtr = eigen_sparse_matrix.innerIndexPtr();
const _Index* outerPtr = eigen_sparse_matrix.outerIndexPtr();
const _Index entryCount = eigen_sparse_matrix.nonZeros();
const _Index outerCount = eigen_sparse_matrix.outerSize();
uint_fast64_t diagNonZeros = 0;
// For RowMajor, row is the current row and col the column and for
// ColMajor, row is the current column and col the row, but this is
// not important as we are only looking for elements on the diagonal.
_Index innerStart = 0;
_Index innerEnd = 0;
_Index innerMid = 0;
for (_Index row = 0; row < outerCount; ++row) {
innerStart = outerPtr[row];
innerEnd = outerPtr[row + 1] - 1;
// Now use binary search (but defer equality detection).
while (innerStart < innerEnd) {
innerMid = innerStart + ((innerEnd - innerStart) / 2);
if (indexPtr[innerMid] < row) {
innerStart = innerMid + 1;
} else {
innerEnd = innerMid;
}
}
// Check whether we have found an element on the diagonal.
if ((innerStart == innerEnd) && (indexPtr[innerStart] == row)) {
++diagNonZeros;
}
}
return static_cast<_Index>(entryCount - diagNonZeros);
}
};
} // namespace storage

113
src/utility/IoUtility.cpp
View File

@ -1,9 +1,9 @@
/*
* IoUtility.cpp
*
* Created on: 17.10.2012
* Author: Thomas Heinemann
*/
* IoUtility.cpp
*
* Created on: 17.10.2012
* Author: Thomas Heinemann
*/
#include "src/utility/IoUtility.h"
#include "src/parser/DeterministicSparseTransitionParser.h"
@ -13,63 +13,58 @@
namespace storm {
namespace utility {
void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename) {
std::shared_ptr<storm::storage::SquareSparseMatrix<double>> matrix(dtmc.getTransitionProbabilityMatrix());
double* diagonal_storage = matrix->getDiagonalStoragePointer();
std::ofstream file;
file.open(filename);
file << "digraph dtmc {\n";
//Specify the nodes and their labels
for (uint_fast64_t i = 1; i < dtmc.getNumberOfStates(); i++) {
file << "\t" << i << "[label=\"" << i << "\\n{";
char komma=' ';
std::set<std::string> propositions = dtmc.getPropositionsForState(i);
for(auto it = propositions.begin();
it != propositions.end();
it++) {
file << komma << *it;
komma=',';
}
file << " }\"];\n";
}
for (uint_fast64_t row = 0; row < dtmc.getNumberOfStates(); row++ ) {
//write diagonal entry/self loop first
if (diagonal_storage[row] != 0) {
file << "\t" << row << " -> " << row << " [label=" << diagonal_storage[row] <<"]\n";
}
//Then, iterate through the row and write each non-diagonal value into the file
for ( auto it = matrix->beginConstColumnNoDiagIterator(row);
it != matrix->endConstColumnNoDiagIterator(row);
it++) {
double value = 0;
matrix->getValue(row,*it,&value);
file << "\t" << row << " -> " << *it << " [label=" << value << "]\n";
}
}
file << "}\n";
file.close();
}
namespace utility {
//TODO: Should this stay here or be integrated in the new parser structure?
/*storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file) {
storm::parser::DeterministicSparseTransitionParser tp(tra_file);
uint_fast64_t node_count = tp.getMatrix()->getRowCount();
void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename) {
std::shared_ptr<storm::storage::SquareSparseMatrix<double>> matrix(dtmc.getTransitionProbabilityMatrix());
std::ofstream file;
file.open(filename);
storm::parser::AtomicPropositionLabelingParser lp(node_count, lab_file);
file << "digraph dtmc {\n";
storm::models::Dtmc<double>* result = new storm::models::Dtmc<double>(tp.getMatrix(), lp.getLabeling());
return result;
}*/
//Specify the nodes and their labels
for (uint_fast64_t i = 1; i < dtmc.getNumberOfStates(); i++) {
file << "\t" << i << "[label=\"" << i << "\\n{";
char komma=' ';
std::set<std::string> propositions = dtmc.getPropositionsForState(i);
for(auto it = propositions.begin();
it != propositions.end();
it++) {
file << komma << *it;
komma=',';
}
}
file << " }\"];\n";
}
for (uint_fast64_t row = 0; row < dtmc.getNumberOfStates(); row++ ) {
//Then, iterate through the row and write each non-diagonal value into the file
for ( auto it = matrix->beginConstColumnIterator(row);
it != matrix->endConstColumnIterator(row);
it++) {
double value = 0;
matrix->getValue(row,*it,&value);
file << "\t" << row << " -> " << *it << " [label=" << value << "]\n";
}
}
file << "}\n";
file.close();
}
//TODO: Should this stay here or be integrated in the new parser structure?
/*storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file) {
storm::parser::DeterministicSparseTransitionParser tp(tra_file);
uint_fast64_t node_count = tp.getMatrix()->getRowCount();
storm::parser::AtomicPropositionLabelingParser lp(node_count, lab_file);
storm::models::Dtmc<double>* result = new storm::models::Dtmc<double>(tp.getMatrix(), lp.getLabeling());
return result;
}*/
}
}

4
test/parser/ReadTraFileTest.cpp
View File

@ -53,13 +53,13 @@ TEST(ReadTraFileTest, ParseFileTest1) {
ASSERT_TRUE(result->getValue(3,2,&val));
ASSERT_EQ(val,0.0806451612903225806451612903225812);
ASSERT_TRUE(result->getValue(3,3,&val));
ASSERT_FALSE(result->getValue(3,3,&val));
ASSERT_EQ(val,0);
ASSERT_TRUE(result->getValue(3,4,&val));
ASSERT_EQ(val,0.080645161290322580645161290322581);
ASSERT_TRUE(result->getValue(4,4,&val));
ASSERT_FALSE(result->getValue(4,4,&val));
ASSERT_EQ(val,0);
delete parser;

20
test/storage/SquareSparseMatrixTest.cpp
View File

@ -116,11 +116,7 @@ TEST(SquareSparseMatrixTest, Test) {
for (int row = 15; row < 24; ++row) {
for (int col = 1; col <= 25; ++col) {
target = 1;
if (row != col) {
ASSERT_FALSE(ssm->getValue(row, col, &target));
} else {
ASSERT_TRUE(ssm->getValue(row, col, &target));
}
ASSERT_FALSE(ssm->getValue(row, col, &target));
ASSERT_EQ(0, target);
}
@ -245,7 +241,7 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest
TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest) {
// 10 rows, 15 non zero entries
storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10, 10);
ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
Eigen::SparseMatrix<int, Eigen::RowMajor> esm(10, 10);
@ -253,7 +249,7 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest
typedef Eigen::Triplet<int> IntTriplet;
std::vector<IntTriplet> tripletList;
tripletList.reserve(15);
tripletList.push_back(IntTriplet(1, 0, 0));
tripletList.push_back(IntTriplet(1, 0, 15));
tripletList.push_back(IntTriplet(1, 1, 1));
tripletList.push_back(IntTriplet(1, 2, 2));
tripletList.push_back(IntTriplet(1, 3, 3));
@ -281,11 +277,15 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest
ASSERT_NO_THROW(ssm->finalize());
ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
const std::vector<uint_fast64_t> rowP = ssm->getRowIndicationsPointer();
const std::vector<uint_fast64_t> colP = ssm->getColumnIndicationsPointer();
const std::vector<int> valP = ssm->getStoragePointer();
int target = -1;
for (auto &coeff: tripletList) {
ASSERT_TRUE(ssm->getValue(coeff.row(), coeff.col(), &target));
bool retVal = ssm->getValue(coeff.row(), coeff.col(), &target);
ASSERT_TRUE(retVal);
ASSERT_EQ(target, coeff.value());
}
@ -300,7 +300,7 @@ TEST(SquareSparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest)
values[i] = i;
}
ASSERT_NO_THROW(ssm->initialize(100 - 10));
ASSERT_NO_THROW(ssm->initialize(100));
ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
for (uint_fast32_t row = 0; row < 10; ++row) {

Write Preview
Loading…
Cancel
Save