From 7bebfb91a01d89c53d68e2cf8e8b59ab84f64f3d Mon Sep 17 00:00:00 2001 From: Stefan Pranger Date: Fri, 29 Jan 2021 08:26:34 +0100 Subject: [PATCH] smg lra debug commit this should be dropped in the future --- src/storm/builder/ExplicitModelBuilder.cpp | 1 + ...deterministicGameInfiniteHorizonHelper.cpp | 3 + .../infinitehorizon/internal/LraViHelper.cpp | 71 ++- src/storm/solver/GmmxxMultiplier.cpp | 77 ++- src/storm/storage/Decomposition.cpp | 56 +- src/storm/storage/MaximalEndComponent.cpp | 66 +- src/storm/storage/SparseMatrix.cpp | 593 +++++++++--------- 7 files changed, 443 insertions(+), 424 deletions(-) diff --git a/src/storm/builder/ExplicitModelBuilder.cpp b/src/storm/builder/ExplicitModelBuilder.cpp index 70b68fbce..3df2f42b7 100644 --- a/src/storm/builder/ExplicitModelBuilder.cpp +++ b/src/storm/builder/ExplicitModelBuilder.cpp @@ -154,6 +154,7 @@ namespace storm { while (!statesToExplore.empty()) { // Get the first state in the queue. CompressedState currentState = statesToExplore.front().first; + STORM_LOG_DEBUG("Exploring (" << currentRowGroup << ") : " << toString(currentState, this->generator->getVariableInformation())); StateType currentIndex = statesToExplore.front().second; statesToExplore.pop_front(); diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicGameInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicGameInfiniteHorizonHelper.cpp index d4f52f65c..34a1a8d93 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicGameInfiniteHorizonHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseNondeterministicGameInfiniteHorizonHelper.cpp @@ -63,7 +63,10 @@ namespace storm { this->_backwardTransitions = this->_computedBackwardTransitions.get(); } this->_computedLongRunComponentDecomposition = std::make_unique>(this->_transitionMatrix, *this->_backwardTransitions); + this->_longRunComponentDecomposition = this->_computedLongRunComponentDecomposition.get(); + //STORM_LOG_DEBUG("\n" << this->_transitionMatrix); + STORM_LOG_DEBUG("GMEC: " << *(this->_longRunComponentDecomposition)); } } diff --git a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp index 5e4bfc81a..da16522f1 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.cpp @@ -24,11 +24,11 @@ namespace storm { namespace modelchecker { namespace helper { namespace internal { - + template LraViHelper::LraViHelper(ComponentType const& component, storm::storage::SparseMatrix const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates, std::vector const* exitRates) : _transitionMatrix(transitionMatrix), _timedStates(timedStates), _hasInstantStates(TransitionsType == LraViTransitionsType::DetTsNondetIs || TransitionsType == LraViTransitionsType::DetTsDetIs), _Tsx1IsCurrent(false) { setComponent(component); - + // Run through the component and collect some data: // We create two submodels, one consisting of the timed states of the component and one consisting of the instant states of the component. // For this, we create a state index map that point from state indices of the input model to indices of the corresponding submodel of that state. @@ -141,8 +141,9 @@ namespace storm { _IsTransitions = isTransitionsBuilder.build(); _IsToTsTransitions = isToTsTransitionsBuilder.build(); } + STORM_LOG_DEBUG(uniformizationFactor << " - " << _uniformizationRate); } - + template void LraViHelper::setComponent(ComponentType component) { _component.clear(); @@ -156,7 +157,7 @@ namespace storm { } } - + template ValueType LraViHelper::performValueIteration(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates, storm::solver::OptimizationDirection const* dir, std::vector* choices) { initializeNewValues(stateValueGetter, actionValueGetter, exitRates); @@ -166,7 +167,7 @@ namespace storm { if (env.solver().lra().isMaximalIterationCountSet()) { maxIter = env.solver().lra().getMaximalIterationCount(); } - + // start the iterations ValueType result = storm::utility::zero(); uint64_t iter = 0; @@ -198,7 +199,7 @@ namespace storm { // If there will be a next iteration, we have to prepare it. prepareNextIteration(env); - + } if (maxIter.is_initialized() && iter == maxIter.get()) { STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations."); @@ -207,15 +208,23 @@ namespace storm { } else { STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations."); } - + if (choices) { // We will be doing one more iteration step and track scheduler choices this time. prepareNextIteration(env); performIterationStep(env, dir, choices); } + std::cout << "result (" << iter << " steps):" << std::endl; + for(int i = 0; i < xNew().size() ; i++ ) { + std::cout << std::setprecision(4) << i << "\t: " << xNew().at(i) << "\t" << xNew().at(i) * _uniformizationRate << "\t" << std::setprecision(16) << xOld().at(i) *_uniformizationRate << std::endl; + //if(i == 50) {std::cout << "only showing top 50 lines"; break; } + for(int i = 0; i < xNew().size() ; i++ ) { + std::cout << std::setprecision(4) << i << "\t: " << xNew().at(i) << "\t" << xNew().at(i) * _uniformizationRate << "\t" << std::setprecision(16) << xOld().at(i) *_uniformizationRate << std::endl; + //if(i == 50) {std::cout << "only showing top 50 lines"; break; } + } return result; } - + template void LraViHelper::initializeNewValues(ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector const* exitRates) { // clear potential old values and reserve enough space for new values @@ -225,7 +234,7 @@ namespace storm { _IsChoiceValues.clear(); _IsChoiceValues.reserve(_IsTransitions.getRowCount()); } - + // Set the new choice-based values ValueType actionRewardScalingFactor = storm::utility::one() / _uniformizationRate; for (auto const& element : _component) { @@ -250,14 +259,14 @@ namespace storm { // Set-up new iteration vectors for timed states _Tsx1.assign(_TsTransitions.getRowGroupCount(), storm::utility::zero()); _Tsx2 = _Tsx1; - + if (_hasInstantStates) { // Set-up vectors for storing intermediate results for instant states. _Isx.resize(_IsTransitions.getRowGroupCount(), storm::utility::zero()); _Isb = _IsChoiceValues; } } - + template void LraViHelper::prepareSolversAndMultipliers(const Environment& env, storm::solver::OptimizationDirection const* dir) { _TsMultiplier = storm::solver::MultiplierFactory().create(env, _TsTransitions); @@ -315,7 +324,7 @@ namespace storm { STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "The solver requirement " << req.getEnabledRequirementsAsString() << " has not been cleared."); } } - + // Set up multipliers for transitions connecting timed and instant states _TsToIsMultiplier = storm::solver::MultiplierFactory().create(env, _TsToIsTransitions); _IsToTsMultiplier = storm::solver::MultiplierFactory().create(env, _IsToTsTransitions); @@ -324,7 +333,7 @@ namespace storm { _TsMultiplier->setOptimizationDirectionOverride(env.solver().multiplier().getOptimizationDirectionOverride().get()); } } - + template void LraViHelper::setInputModelChoices(std::vector& choices, std::vector const& localMecChoices, bool setChoiceZeroToTimedStates, bool setChoiceZeroToInstantStates) const { // Transform the local choices (within this mec) to choice indices for the input model @@ -347,7 +356,7 @@ namespace storm { } STORM_LOG_ASSERT(localState == localMecChoices.size(), "Did not traverse all component states."); } - + template void LraViHelper::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const* dir, std::vector* choices) { STORM_LOG_ASSERT(!((nondetTs() || nondetIs()) && dir == nullptr), "No optimization direction provided for model with nondeterminism"); @@ -355,13 +364,13 @@ namespace storm { if (!_TsMultiplier) { prepareSolversAndMultipliers(env, dir); } - + // Compute new x values for the timed states // Flip what is new and what is old _Tsx1IsCurrent = !_Tsx1IsCurrent; // At this point, xOld() points to what has been computed in the most recent call of performIterationStep (initially, this is the 0-vector). // The result of this ongoing computation will be stored in xNew() - + // Compute the values obtained by a single uniformization step between timed states only if (nondetTs() && !gameNondetTs()) { if (choices == nullptr) { @@ -425,7 +434,7 @@ namespace storm { _TsToIsMultiplier->multiply(env, _Isx, &xNew(), xNew()); } } - + template typename LraViHelper::ConvergenceCheckResult LraViHelper::checkConvergence(bool relative, ValueType precision) const { STORM_LOG_ASSERT(_TsMultiplier, "tried to check for convergence without doing an iteration first."); @@ -433,7 +442,7 @@ namespace storm { // We need to 'revert' this scaling when computing the absolute precision. // However, for relative precision, the scaling cancels out. ValueType threshold = relative ? precision : ValueType(precision / _uniformizationRate); - + ConvergenceCheckResult res = { true, storm::utility::one() }; // Now check whether the currently produced results are precise enough STORM_LOG_ASSERT(threshold > storm::utility::zero(), "Did not expect a non-positive threshold."); @@ -460,15 +469,15 @@ namespace storm { break; } } - + // Compute the average of the maximal and the minimal difference. ValueType avgDiff = (maxDiff + minDiff) / (storm::utility::convertNumber(2.0)); - + // "Undo" the scaling of the values res.currentValue = avgDiff * _uniformizationRate; return res; } - + template void LraViHelper::prepareNextIteration(Environment const& env) { // To avoid large (and numerically unstable) x-values, we substract a reference value. @@ -480,39 +489,39 @@ namespace storm { _IsToTsMultiplier->multiply(env, xNew(), &_IsChoiceValues, _Isb); } } - + template bool LraViHelper::isTimedState(uint64_t const& inputModelStateIndex) const { STORM_LOG_ASSERT(!_hasInstantStates || _timedStates != nullptr, "Model has instant states but no partition into timed and instant states is given."); STORM_LOG_ASSERT(!_hasInstantStates || inputModelStateIndex < _timedStates->size(), "Unable to determine whether state " << inputModelStateIndex << " is timed."); return !_hasInstantStates || _timedStates->get(inputModelStateIndex); } - + template std::vector& LraViHelper::xNew() { return _Tsx1IsCurrent ? _Tsx1 : _Tsx2; } - + template std::vector const& LraViHelper::xNew() const { return _Tsx1IsCurrent ? _Tsx1 : _Tsx2; } - + template std::vector& LraViHelper::xOld() { return _Tsx1IsCurrent ? _Tsx2 : _Tsx1; } - + template std::vector const& LraViHelper::xOld() const { return _Tsx1IsCurrent ? _Tsx2 : _Tsx1; } - + template bool LraViHelper::nondetTs() const { return TransitionsType == LraViTransitionsType::NondetTsNoIs || gameNondetTs(); } - + template bool LraViHelper::nondetIs() const { return TransitionsType == LraViTransitionsType::DetTsNondetIs; @@ -529,11 +538,11 @@ namespace storm { template class LraViHelper; template class LraViHelper; template class LraViHelper; - + template class LraViHelper; template class LraViHelper; - + } } } -} \ No newline at end of file +} diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp index 32af42f72..0a2843824 100644 --- a/src/storm/solver/GmmxxMultiplier.cpp +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -16,25 +16,26 @@ namespace storm { namespace solver { - + template GmmxxMultiplier::GmmxxMultiplier(storm::storage::SparseMatrix const& matrix) : Multiplier(matrix) { // Intentionally left empty. + //STORM_LOG_DEBUG("\n" << matrix); } - + template void GmmxxMultiplier::initialize() const { if (gmmMatrix.nrows() == 0) { gmmMatrix = std::move(*storm::adapters::GmmxxAdapter().toGmmxxSparseMatrix(this->matrix)); } } - + template void GmmxxMultiplier::clearCache() const { gmmMatrix = gmm::csr_matrix(); Multiplier::clearCache(); } - + template bool GmmxxMultiplier::parallelize(Environment const& env) const { #ifdef STORM_HAVE_INTELTBB @@ -43,7 +44,7 @@ namespace storm { return false; #endif } - + template void GmmxxMultiplier::multiply(Environment const& env, std::vector const& x, std::vector const* b, std::vector& result) const { initialize(); @@ -65,7 +66,7 @@ namespace storm { std::swap(result, *this->cachedVector); } } - + template void GmmxxMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards) const { initialize(); @@ -84,7 +85,7 @@ namespace storm { } } } - + template void GmmxxMultiplier::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { initialize(); @@ -106,13 +107,13 @@ namespace storm { std::swap(result, *this->cachedVector); } } - + template void GmmxxMultiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices, bool backwards) const { initialize(); multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices, backwards); } - + template void GmmxxMultiplier::multiplyRow(uint64_t const& rowIndex, std::vector const& x, ValueType& value) const { initialize(); @@ -148,14 +149,14 @@ namespace storm { } } } - + template template void GmmxxMultiplier::multAddReduceHelper(std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { Compare compare; typedef std::vector VectorType; typedef gmm::csr_matrix MatrixType; - + typename gmm::linalg_traits::const_iterator add_it, add_ite; if (b) { add_it = backwards ? gmm::vect_end(*b) - 1 : gmm::vect_begin(*b); @@ -181,6 +182,7 @@ namespace storm { uint64_t currentRowGroup = backwards ? rowGroupIndices.size() - 1 : 0; auto row_group_it = backwards ? rowGroupIndices.end() - 2 : rowGroupIndices.begin(); auto row_group_ite = backwards ? rowGroupIndices.begin() - 1 : rowGroupIndices.end() - 1; + //if(choices) STORM_LOG_DEBUG(" "); while (row_group_it != row_group_ite) { ValueType currentValue = storm::utility::zero(); @@ -199,7 +201,7 @@ namespace storm { oldSelectedChoiceValue = currentValue; } } - + // move row-based iterators to the next row if (backwards) { --itr; @@ -213,6 +215,10 @@ namespace storm { // Process the (rowGroupSize-1) remaining rows within the current row Group uint64_t rowGroupSize = *(row_group_it + 1) - *row_group_it; + uint choiceforprintout = 0; + //std::cout << currentRowGroup << ": " << currentValue << ", "; + //STORM_LOG_DEBUG(std::setprecision(3) << vect_sp(gmm::linalg_traits::row(itr), x) << " + " << *add_it << "; "); + //STORM_LOG_DEBUG(std::setprecision(3) << vect_sp(gmm::linalg_traits::row(itr), x) << " + " << *add_it << "; "); for (uint64_t i = 1; i < rowGroupSize; ++i) { ValueType newValue = b ? *add_it : storm::utility::zero(); newValue += vect_sp(gmm::linalg_traits::row(itr), x); @@ -220,12 +226,13 @@ namespace storm { if (choices && currentRow == *choice_it + *row_group_it) { oldSelectedChoiceValue = newValue; } - - if(this->isOverridden(currentRowGroup) ? !compare(newValue, currentValue) : compare(newValue, currentValue)) { + //std::cout << newValue << ", "; + //STORM_LOG_DEBUG(std::setprecision(3) << vect_sp(gmm::linalg_traits::row(itr), x) << " + " << *add_it << "; "); currentValue = newValue; if (choices) { selectedChoice = currentRow - *row_group_it; } + choiceforprintout = currentRow - *row_group_it; } // move row-based iterators to the next row if (backwards) { @@ -238,7 +245,8 @@ namespace storm { ++add_it; } } - + //STORM_LOG_DEBUG("\t= " << currentValue << "\tchoice: " << choiceforprintout); + //std::cout << std::fixed << std::setprecision(2) << " | v(" << currentRowGroup << ")=" << currentValue << " c: " << choiceforprintout << " |\n" ; // Finally write value to target vector. *target_it = currentValue; if(choices) { @@ -261,6 +269,7 @@ namespace storm { ++currentRowGroup; } } + //std::cout << std::endl; } template<> @@ -286,7 +295,7 @@ namespace storm { multAdd(x, b, result); #endif } - + #ifdef STORM_HAVE_INTELTBB template class TbbMultAddReduceFunctor { @@ -294,14 +303,14 @@ namespace storm { TbbMultAddReduceFunctor(std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) : rowGroupIndices(rowGroupIndices), matrix(matrix), x(x), b(b), result(result), choices(choices) { // Intentionally left empty. } - + void operator()(tbb::blocked_range const& range) const { typedef std::vector VectorType; typedef gmm::csr_matrix MatrixType; auto groupIt = rowGroupIndices.begin() + range.begin(); auto groupIte = rowGroupIndices.begin() + range.end(); - + auto itr = mat_row_const_begin(matrix) + *groupIt; typename std::vector::const_iterator bIt; if (b) { @@ -311,40 +320,40 @@ namespace storm { if (choices) { choiceIt = choices->begin() + range.begin(); } - + auto resultIt = result.begin() + range.begin(); - + // Variables for correctly tracking choices (only update if new choice is strictly better). ValueType oldSelectedChoiceValue; uint64_t selectedChoice; - + uint64_t currentRow = *groupIt; for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { ValueType currentValue = storm::utility::zero(); - + // Only multiply and reduce if the row group is not empty. if (*groupIt != *(groupIt + 1)) { if (b) { currentValue = *bIt; ++bIt; } - + currentValue += vect_sp(gmm::linalg_traits::row(itr), x); - + if (choices) { selectedChoice = currentRow - *groupIt; if (*choiceIt == selectedChoice) { oldSelectedChoiceValue = currentValue; } } - + ++itr; ++currentRow; - + for (auto itre = mat_row_const_begin(matrix) + *(groupIt + 1); itr != itre; ++itr, ++bIt, ++currentRow) { ValueType newValue = b ? *bIt : storm::utility::zero(); newValue += vect_sp(gmm::linalg_traits::row(itr), x); - + if (compare(newValue, currentValue)) { currentValue = newValue; if (choices) { @@ -353,7 +362,7 @@ namespace storm { } } } - + // Finally write value to target vector. *resultIt = currentValue; if (choices && compare(currentValue, oldSelectedChoiceValue)) { @@ -361,7 +370,7 @@ namespace storm { } } } - + private: Compare compare; std::vector const& rowGroupIndices; @@ -372,7 +381,7 @@ namespace storm { std::vector* choices; }; #endif - + template void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { #ifdef STORM_HAVE_INTELTBB @@ -386,18 +395,18 @@ namespace storm { multAddReduceHelper(dir, rowGroupIndices, x, b, result, choices); #endif } - + template<> void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } - + template class GmmxxMultiplier; - + #ifdef STORM_HAVE_CARL template class GmmxxMultiplier; template class GmmxxMultiplier; #endif - + } } diff --git a/src/storm/storage/Decomposition.cpp b/src/storm/storage/Decomposition.cpp index efe850083..a45987d69 100644 --- a/src/storm/storage/Decomposition.cpp +++ b/src/storm/storage/Decomposition.cpp @@ -8,84 +8,84 @@ namespace storm { namespace storage { - + template Decomposition::Decomposition() : blocks() { // Intentionally left empty. } - + template Decomposition::Decomposition(Decomposition const& other) : blocks(other.blocks) { // Intentionally left empty. } - + template Decomposition& Decomposition::operator=(Decomposition const& other) { this->blocks = other.blocks; return *this; } - + template Decomposition::Decomposition(Decomposition&& other) : blocks(std::move(other.blocks)) { // Intentionally left empty. } - + template Decomposition& Decomposition::operator=(Decomposition&& other) { this->blocks = std::move(other.blocks); return *this; } - + template std::size_t Decomposition::size() const { return blocks.size(); } - + template bool Decomposition::empty() const { return blocks.empty(); } - + template typename Decomposition::iterator Decomposition::begin() { return blocks.begin(); } - + template typename Decomposition::iterator Decomposition::end() { return blocks.end(); } - + template typename Decomposition::const_iterator Decomposition::begin() const { return blocks.begin(); } - + template typename Decomposition::const_iterator Decomposition::end() const { return blocks.end(); } - + template BlockType const& Decomposition::getBlock(uint_fast64_t index) const { return blocks.at(index); } - + template BlockType& Decomposition::getBlock(uint_fast64_t index) { return blocks.at(index); } - + template BlockType const& Decomposition::operator[](uint_fast64_t index) const { return blocks[index]; } - + template BlockType& Decomposition::operator[](uint_fast64_t index) { return blocks[index]; } - + template template storm::storage::SparseMatrix Decomposition::extractPartitionDependencyGraph(storm::storage::SparseMatrix const& matrix) const { @@ -97,49 +97,49 @@ namespace storm { stateToBlockMap[state] = i; } } - + // The resulting sparse matrix will have as many rows/columns as there are blocks in the partition. storm::storage::SparseMatrixBuilder dependencyGraphBuilder(this->size(), this->size()); - + for (uint_fast64_t currentBlockIndex = 0; currentBlockIndex < this->size(); ++currentBlockIndex) { // Get the next block. block_type const& block = this->getBlock(currentBlockIndex); - + // Now, we determine the blocks which are reachable (in one step) from the current block. storm::storage::FlatSet allTargetBlocks; for (auto state : block) { for (auto const& transitionEntry : matrix.getRowGroup(state)) { uint_fast64_t targetBlock = stateToBlockMap[transitionEntry.getColumn()]; - + // We only need to consider transitions that are actually leaving the SCC. if (targetBlock != currentBlockIndex) { allTargetBlocks.insert(targetBlock); } } } - + // Now we can just enumerate all the target blocks and insert the corresponding transitions. for (auto const& targetBlock : allTargetBlocks) { dependencyGraphBuilder.addNextValue(currentBlockIndex, targetBlock, storm::utility::one()); } } - + return dependencyGraphBuilder.build(); } - + template std::ostream& operator<<(std::ostream& out, Decomposition const& decomposition) { - out << "["; + out << "[ "; if (decomposition.size() > 0) { for (uint_fast64_t blockIndex = 0; blockIndex < decomposition.size() - 1; ++blockIndex) { - out << decomposition.blocks[blockIndex] << ", "; + out << decomposition.blocks[blockIndex] << ", " << std::endl; } out << decomposition.blocks.back(); } - out << "]"; + out << " ]"; return out; } - + template storm::storage::SparseMatrix Decomposition::extractPartitionDependencyGraph(storm::storage::SparseMatrix const& matrix) const; template storm::storage::SparseMatrix Decomposition::extractPartitionDependencyGraph(storm::storage::SparseMatrix const& matrix) const; template class Decomposition; @@ -149,7 +149,7 @@ namespace storm { template storm::storage::SparseMatrix Decomposition::extractPartitionDependencyGraph(storm::storage::SparseMatrix const& matrix) const; template class Decomposition; template std::ostream& operator<<(std::ostream& out, Decomposition const& decomposition); - + template class Decomposition; template std::ostream& operator<<(std::ostream& out, Decomposition const& decomposition); } // namespace storage diff --git a/src/storm/storage/MaximalEndComponent.cpp b/src/storm/storage/MaximalEndComponent.cpp index 68e8eda99..58683a3df 100644 --- a/src/storm/storage/MaximalEndComponent.cpp +++ b/src/storm/storage/MaximalEndComponent.cpp @@ -3,125 +3,125 @@ namespace storm { namespace storage { - + std::ostream& operator<<(std::ostream& out, storm::storage::FlatSet const& block); - + MaximalEndComponent::MaximalEndComponent() : stateToChoicesMapping() { // Intentionally left empty. } - + MaximalEndComponent::MaximalEndComponent(MaximalEndComponent const& other) : stateToChoicesMapping(other.stateToChoicesMapping) { // Intentionally left empty. } - + MaximalEndComponent& MaximalEndComponent::operator=(MaximalEndComponent const& other) { stateToChoicesMapping = other.stateToChoicesMapping; return *this; } - + MaximalEndComponent::MaximalEndComponent(MaximalEndComponent&& other) : stateToChoicesMapping(std::move(other.stateToChoicesMapping)) { // Intentionally left empty. } - + MaximalEndComponent& MaximalEndComponent::operator=(MaximalEndComponent&& other) { stateToChoicesMapping = std::move(other.stateToChoicesMapping); return *this; } - + void MaximalEndComponent::addState(uint_fast64_t state, set_type const& choices) { stateToChoicesMapping[state] = choices; } - + void MaximalEndComponent::addState(uint_fast64_t state, set_type&& choices) { stateToChoicesMapping.emplace(state, std::move(choices)); } - + std::size_t MaximalEndComponent::size() const { return stateToChoicesMapping.size(); } - + MaximalEndComponent::set_type const& MaximalEndComponent::getChoicesForState(uint_fast64_t state) const { auto stateChoicePair = stateToChoicesMapping.find(state); - + if (stateChoicePair == stateToChoicesMapping.end()) { throw storm::exceptions::InvalidStateException() << "Invalid call to MaximalEndComponent::getChoicesForState: cannot retrieve choices for state not contained in MEC."; } - + return stateChoicePair->second; } - + MaximalEndComponent::set_type& MaximalEndComponent::getChoicesForState(uint_fast64_t state) { auto stateChoicePair = stateToChoicesMapping.find(state); - + if (stateChoicePair == stateToChoicesMapping.end()) { throw storm::exceptions::InvalidStateException() << "Invalid call to MaximalEndComponent::getChoicesForState: cannot retrieve choices for state not contained in MEC."; } - + return stateChoicePair->second; } - + bool MaximalEndComponent::containsState(uint_fast64_t state) const { auto stateChoicePair = stateToChoicesMapping.find(state); - + if (stateChoicePair == stateToChoicesMapping.end()) { return false; } return true; } - + void MaximalEndComponent::removeState(uint_fast64_t state) { auto stateChoicePair = stateToChoicesMapping.find(state); - + if (stateChoicePair == stateToChoicesMapping.end()) { throw storm::exceptions::InvalidStateException() << "Invalid call to MaximalEndComponent::removeState: cannot remove state not contained in MEC."; } - + stateToChoicesMapping.erase(stateChoicePair); } - + bool MaximalEndComponent::containsChoice(uint_fast64_t state, uint_fast64_t choice) const { auto stateChoicePair = stateToChoicesMapping.find(state); - + if (stateChoicePair == stateToChoicesMapping.end()) { throw storm::exceptions::InvalidStateException() << "Invalid call to MaximalEndComponent::containsChoice: cannot obtain choices for state not contained in MEC."; } - + return stateChoicePair->second.find(choice) != stateChoicePair->second.end(); } - + MaximalEndComponent::set_type MaximalEndComponent::getStateSet() const { set_type states; states.reserve(stateToChoicesMapping.size()); - + for (auto const& stateChoicesPair : stateToChoicesMapping) { states.insert(stateChoicesPair.first); } - + return states; } - + std::ostream& operator<<(std::ostream& out, MaximalEndComponent const& component) { out << "{"; for (auto const& stateChoicesPair : component.stateToChoicesMapping) { - out << "{" << stateChoicesPair.first << ", " << stateChoicesPair.second << "}"; + out << "(" << stateChoicesPair.first << ", " << stateChoicesPair.second << ")"; } out << "}"; - + return out; } - + MaximalEndComponent::iterator MaximalEndComponent::begin() { return stateToChoicesMapping.begin(); } - + MaximalEndComponent::iterator MaximalEndComponent::end() { return stateToChoicesMapping.end(); } - + MaximalEndComponent::const_iterator MaximalEndComponent::begin() const { return stateToChoicesMapping.begin(); } - + MaximalEndComponent::const_iterator MaximalEndComponent::end() const { return stateToChoicesMapping.end(); } diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 05a526702..d89f7a568 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -21,64 +21,64 @@ namespace storm { namespace storage { - + template MatrixEntry::MatrixEntry(IndexType column, ValueType value) : entry(column, value) { // Intentionally left empty. } - + template MatrixEntry::MatrixEntry(std::pair&& pair) : entry(std::move(pair)) { // Intentionally left empty. } - + template IndexType const& MatrixEntry::getColumn() const { return this->entry.first; } - + template void MatrixEntry::setColumn(IndexType const& column) { this->entry.first = column; } - + template ValueType const& MatrixEntry::getValue() const { return this->entry.second; } - + template void MatrixEntry::setValue(ValueType const& value) { this->entry.second = value; } - + template std::pair const& MatrixEntry::getColumnValuePair() const { return this->entry; } - + template MatrixEntry MatrixEntry::operator*(value_type factor) const { return MatrixEntry(this->getColumn(), this->getValue() * factor); } - - + + template bool MatrixEntry::operator==(MatrixEntry const& other) const { return this->entry.first == other.entry.first && this->entry.second == other.entry.second; } - + template bool MatrixEntry::operator!=(MatrixEntry const& other) const { return !(*this == other); } - + template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry) { out << "(" << entry.getColumn() << ", " << entry.getValue() << ")"; return out; } - + template SparseMatrixBuilder::SparseMatrixBuilder(index_type rows, index_type columns, index_type entries, bool forceDimensions, bool hasCustomRowGrouping, index_type rowGroups) : initialRowCountSet(rows != 0), initialRowCount(rows), initialColumnCountSet(columns != 0), initialColumnCount(columns), initialEntryCountSet(entries != 0), initialEntryCount(entries), forceInitialDimensions(forceDimensions), hasCustomRowGrouping(hasCustomRowGrouping), initialRowGroupCountSet(rowGroups != 0), initialRowGroupCount(rowGroups), rowGroupIndices(), columnsAndValues(), rowIndications(), currentEntryCount(0), lastRow(0), lastColumn(0), highestColumn(0), currentRowGroupCount(0) { // Prepare the internal storage. @@ -96,14 +96,14 @@ namespace storm { } rowIndications.push_back(0); } - + template SparseMatrixBuilder::SparseMatrixBuilder(SparseMatrix&& matrix) : initialRowCountSet(false), initialRowCount(0), initialColumnCountSet(false), initialColumnCount(0), initialEntryCountSet(false), initialEntryCount(0), forceInitialDimensions(false), hasCustomRowGrouping(!matrix.trivialRowGrouping), initialRowGroupCountSet(false), initialRowGroupCount(0), rowGroupIndices(), columnsAndValues(std::move(matrix.columnsAndValues)), rowIndications(std::move(matrix.rowIndications)), currentEntryCount(matrix.entryCount), currentRowGroupCount() { - + lastRow = matrix.rowCount == 0 ? 0 : matrix.rowCount - 1; lastColumn = columnsAndValues.empty() ? 0 : columnsAndValues.back().getColumn(); highestColumn = matrix.getColumnCount() == 0 ? 0 : matrix.getColumnCount() - 1; - + // If the matrix has a custom row grouping, we move it and remove the last element to make it 'open' again. if (hasCustomRowGrouping) { rowGroupIndices = std::move(matrix.rowGroupIndices); @@ -112,19 +112,19 @@ namespace storm { } currentRowGroupCount = rowGroupIndices->empty() ? 0 : rowGroupIndices.get().size() - 1; } - + // Likewise, we need to 'open' the row indications again. if (!rowIndications.empty()) { rowIndications.pop_back(); } } - + template void SparseMatrixBuilder::addNextValue(index_type row, index_type column, ValueType const& value) { // Check that we did not move backwards wrt. the row. STORM_LOG_THROW(row >= lastRow, storm::exceptions::InvalidArgumentException, "Adding an element in row " << row << ", but an element in row " << lastRow << " has already been added."); STORM_LOG_ASSERT(columnsAndValues.size() == currentEntryCount, "Unexpected size of columnsAndValues vector."); - + // Check if a diagonal entry shall be inserted before if (pendingDiagonalEntry) { index_type diagColumn = hasCustomRowGrouping ? currentRowGroupCount - 1 : lastRow; @@ -143,7 +143,7 @@ namespace storm { } } } - + // If the element is in the same row, but was not inserted in the correct order, we need to fix the row after // the insertion. bool fixCurrentRow = row == lastRow && column < lastColumn; @@ -159,26 +159,26 @@ namespace storm { rowIndications.resize(row + 1, currentEntryCount); lastRow = row; } - + lastColumn = column; - + // Finally, set the element and increase the current size. columnsAndValues.emplace_back(column, value); highestColumn = std::max(highestColumn, column); ++currentEntryCount; - + // If we need to fix the row, do so now. if (fixCurrentRow) { // First, we sort according to columns. std::sort(columnsAndValues.begin() + rowIndications.back(), columnsAndValues.end(), [] (storm::storage::MatrixEntry const& a, storm::storage::MatrixEntry const& b) { return a.getColumn() < b.getColumn(); }); - + // Then, we eliminate possible duplicate entries. auto it = std::unique(columnsAndValues.begin() + rowIndications.back(), columnsAndValues.end(), [] (storm::storage::MatrixEntry const& a, storm::storage::MatrixEntry const& b) { return a.getColumn() == b.getColumn(); }); - + // Finally, remove the superfluous elements. std::size_t elementsToRemove = std::distance(it, columnsAndValues.end()); if (elementsToRemove > 0) { @@ -188,7 +188,7 @@ namespace storm { } } } - + // In case we did not expect this value, we throw an exception. if (forceInitialDimensions) { STORM_LOG_THROW(!initialRowCountSet || lastRow < initialRowCount, storm::exceptions::OutOfRangeException, "Cannot insert value at illegal row " << lastRow << "."); @@ -196,12 +196,12 @@ namespace storm { STORM_LOG_THROW(!initialEntryCountSet || currentEntryCount <= initialEntryCount, storm::exceptions::OutOfRangeException, "Too many entries in matrix, expected only " << initialEntryCount << "."); } } - + template void SparseMatrixBuilder::newRowGroup(index_type startingRow) { STORM_LOG_THROW(hasCustomRowGrouping, storm::exceptions::InvalidStateException, "Matrix was not created to have a custom row grouping."); STORM_LOG_THROW(startingRow >= lastRow, storm::exceptions::InvalidStateException, "Illegal row group with negative size."); - + // If there still is a pending diagonal entry, we need to add it now (otherwise, the correct diagonal column will be unclear) if (pendingDiagonalEntry) { STORM_LOG_ASSERT(currentRowGroupCount > 0, "Diagonal entry was set before opening the first row group."); @@ -210,10 +210,10 @@ namespace storm { pendingDiagonalEntry = boost::none; // clear now, so addNextValue works properly addNextValue(lastRow, diagColumn, diagValue); } - + rowGroupIndices.get().push_back(startingRow); ++currentRowGroupCount; - + // Handle the case where the previous row group ends with one or more empty rows if (lastRow + 1 < startingRow) { // Close all rows from the most recent one to the starting row. @@ -224,10 +224,10 @@ namespace storm { lastColumn = 0; } } - + template SparseMatrix SparseMatrixBuilder::build(index_type overriddenRowCount, index_type overriddenColumnCount, index_type overriddenRowGroupCount) { - + // If there still is a pending diagonal entry, we need to add it now if (pendingDiagonalEntry) { index_type diagColumn = hasCustomRowGrouping ? currentRowGroupCount - 1 : lastRow; @@ -235,30 +235,30 @@ namespace storm { pendingDiagonalEntry = boost::none; // clear now, so addNextValue works properly addNextValue(lastRow, diagColumn, diagValue); } - + bool hasEntries = currentEntryCount != 0; - + uint_fast64_t rowCount = hasEntries ? lastRow + 1 : 0; - + // If the last row group was empty, we need to add one more to the row count, because otherwise this empty row is not counted. if (hasCustomRowGrouping) { if (lastRow < rowGroupIndices->back()) { ++rowCount; } } - + if (initialRowCountSet && forceInitialDimensions) { STORM_LOG_THROW(rowCount <= initialRowCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialRowCount << " rows, but got " << rowCount << "."); rowCount = std::max(rowCount, initialRowCount); } - + rowCount = std::max(rowCount, overriddenRowCount); - + // If the current row count was overridden, we may need to add empty rows. for (index_type i = lastRow + 1; i < rowCount; ++i) { rowIndications.push_back(currentEntryCount); } - + // We put a sentinel element at the last position of the row indices array. This eases iteration work, // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for // the first and last row. @@ -272,12 +272,12 @@ namespace storm { columnCount = std::max(columnCount, initialColumnCount); } columnCount = std::max(columnCount, overriddenColumnCount); - + uint_fast64_t entryCount = currentEntryCount; if (initialEntryCountSet && forceInitialDimensions) { STORM_LOG_THROW(entryCount == initialEntryCount, storm::exceptions::InvalidStateException, "Expected " << initialEntryCount << " entries, but got " << entryCount << "."); } - + // Check whether row groups are missing some entries. if (hasCustomRowGrouping) { uint_fast64_t rowGroupCount = currentRowGroupCount; @@ -286,15 +286,15 @@ namespace storm { rowGroupCount = std::max(rowGroupCount, initialRowGroupCount); } rowGroupCount = std::max(rowGroupCount, overriddenRowGroupCount); - + for (index_type i = currentRowGroupCount; i <= rowGroupCount; ++i) { rowGroupIndices.get().push_back(rowCount); } } - + return SparseMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)); } - + template typename SparseMatrixBuilder::index_type SparseMatrixBuilder::getLastRow() const { return lastRow; @@ -313,7 +313,7 @@ namespace storm { typename SparseMatrixBuilder::index_type SparseMatrixBuilder::getLastColumn() const { return lastColumn; } - + // Debug method for printing the current matrix template void print(std::vector::index_type> const& rowGroupIndices, std::vector::index_type, typename SparseMatrix::value_type>> const& columnsAndValues, std::vector::index_type> const& rowIndications) { @@ -335,11 +335,11 @@ namespace storm { } } } - + template void SparseMatrixBuilder::replaceColumns(std::vector const& replacements, index_type offset) { index_type maxColumn = 0; - + for (index_type row = 0; row < rowIndications.size(); ++row) { bool changed = false; auto startRow = std::next(columnsAndValues.begin(), rowIndications[row]); @@ -365,11 +365,11 @@ namespace storm { }), "Columns not sorted."); } } - + highestColumn = maxColumn; lastColumn = columnsAndValues.empty() ? 0 : columnsAndValues.back().getColumn(); } - + template void SparseMatrixBuilder::addDiagonalEntry(index_type row, ValueType const& value) { STORM_LOG_THROW(row >= lastRow, storm::exceptions::InvalidArgumentException, "Adding a diagonal element in row " << row << ", but an element in row " << lastRow << " has already been added."); @@ -399,59 +399,59 @@ namespace storm { SparseMatrix::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) { // Intentionally left empty. } - + template typename SparseMatrix::iterator SparseMatrix::rows::begin() { return beginIterator; } - + template typename SparseMatrix::iterator SparseMatrix::rows::end() { return beginIterator + entryCount; } - + template typename SparseMatrix::index_type SparseMatrix::rows::getNumberOfEntries() const { return this->entryCount; } - + template SparseMatrix::const_rows::const_rows(const_iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) { // Intentionally left empty. } - + template typename SparseMatrix::const_iterator SparseMatrix::const_rows::begin() const { return beginIterator; } - + template typename SparseMatrix::const_iterator SparseMatrix::const_rows::end() const { return beginIterator + entryCount; } - + template typename SparseMatrix::index_type SparseMatrix::const_rows::getNumberOfEntries() const { return this->entryCount; } - + template SparseMatrix::SparseMatrix() : rowCount(0), columnCount(0), entryCount(0), nonzeroEntryCount(0), columnsAndValues(), rowIndications(), rowGroupIndices() { // Intentionally left empty. } - + template SparseMatrix::SparseMatrix(SparseMatrix const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), trivialRowGrouping(other.trivialRowGrouping), rowGroupIndices(other.rowGroupIndices) { // Intentionally left empty. } - + template SparseMatrix::SparseMatrix(SparseMatrix const& other, bool insertDiagonalElements) { storm::storage::BitVector rowConstraint(other.getRowCount(), true); storm::storage::BitVector columnConstraint(other.getColumnCount(), true); *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements); } - + template SparseMatrix::SparseMatrix(SparseMatrix&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), trivialRowGrouping(other.trivialRowGrouping), rowGroupIndices(std::move(other.rowGroupIndices)) { // Now update the source matrix @@ -459,12 +459,12 @@ namespace storm { other.columnCount = 0; other.entryCount = 0; } - + template SparseMatrix::SparseMatrix(index_type columnCount, std::vector const& rowIndications, std::vector> const& columnsAndValues, boost::optional> const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), trivialRowGrouping(!rowGroupIndices), rowGroupIndices(rowGroupIndices) { this->updateNonzeroEntryCount(); } - + template SparseMatrix::SparseMatrix(index_type columnCount, std::vector&& rowIndications, std::vector>&& columnsAndValues, boost::optional>&& rowGroupIndices) : columnCount(columnCount), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) { // Initialize some variables here which depend on other variables @@ -474,7 +474,7 @@ namespace storm { this->trivialRowGrouping = !this->rowGroupIndices; this->updateNonzeroEntryCount(); } - + template SparseMatrix& SparseMatrix::operator=(SparseMatrix const& other) { // Only perform assignment if source and target are not the same. @@ -483,7 +483,7 @@ namespace storm { columnCount = other.columnCount; entryCount = other.entryCount; nonzeroEntryCount = other.nonzeroEntryCount; - + columnsAndValues = other.columnsAndValues; rowIndications = other.rowIndications; rowGroupIndices = other.rowGroupIndices; @@ -491,7 +491,7 @@ namespace storm { } return *this; } - + template SparseMatrix& SparseMatrix::operator=(SparseMatrix&& other) { // Only perform assignment if source and target are not the same. @@ -500,7 +500,7 @@ namespace storm { columnCount = other.columnCount; entryCount = other.entryCount; nonzeroEntryCount = other.nonzeroEntryCount; - + columnsAndValues = std::move(other.columnsAndValues); rowIndications = std::move(other.rowIndications); rowGroupIndices = std::move(other.rowGroupIndices); @@ -508,15 +508,15 @@ namespace storm { } return *this; } - + template bool SparseMatrix::operator==(SparseMatrix const& other) const { if (this == &other) { return true; } - + bool equalityResult = true; - + equalityResult &= this->getRowCount() == other.getRowCount(); if (!equalityResult) { return false; @@ -533,7 +533,7 @@ namespace storm { if (!equalityResult) { return false; } - + // For the actual contents, we need to do a little bit more work, because we want to ignore elements that // are set to zero, please they may be represented implicitly in the other matrix. for (index_type row = 0; row < this->getRowCount(); ++row) { @@ -559,25 +559,25 @@ namespace storm { return false; } } - + return equalityResult; } - + template typename SparseMatrix::index_type SparseMatrix::getRowCount() const { return rowCount; } - + template typename SparseMatrix::index_type SparseMatrix::getColumnCount() const { return columnCount; } - + template typename SparseMatrix::index_type SparseMatrix::getEntryCount() const { return entryCount; } - + template uint_fast64_t SparseMatrix::getRowGroupEntryCount(uint_fast64_t const group) const { uint_fast64_t result = 0; @@ -590,12 +590,12 @@ namespace storm { } return result; } - + template typename SparseMatrix::index_type SparseMatrix::getNonzeroEntryCount() const { return nonzeroEntryCount; } - + template void SparseMatrix::updateNonzeroEntryCount() const { this->nonzeroEntryCount = 0; @@ -605,12 +605,12 @@ namespace storm { } } } - + template void SparseMatrix::updateNonzeroEntryCount(std::make_signed::type difference) { this->nonzeroEntryCount += difference; } - + template void SparseMatrix::updateDimensions() const { this->nonzeroEntryCount = 0; @@ -622,7 +622,7 @@ namespace storm { } } } - + template typename SparseMatrix::index_type SparseMatrix::getRowGroupCount() const { if (!this->hasTrivialRowGrouping()) { @@ -631,12 +631,12 @@ namespace storm { return rowCount; } } - + template typename SparseMatrix::index_type SparseMatrix::getRowGroupSize(index_type group) const { return this->getRowGroupIndices()[group + 1] - this->getRowGroupIndices()[group]; } - + template typename SparseMatrix::index_type SparseMatrix::getSizeOfLargestRowGroup() const { if (this->hasTrivialRowGrouping()) { @@ -650,7 +650,7 @@ namespace storm { } return res; } - + template typename SparseMatrix::index_type SparseMatrix::getNumRowsInRowGroups(storm::storage::BitVector const& groupConstraint) const { if (this->hasTrivialRowGrouping()) { @@ -669,7 +669,7 @@ namespace storm { return numRows; } - + template std::vector::index_type> const& SparseMatrix::getRowGroupIndices() const { // If there is no current row grouping, we need to create it. @@ -679,7 +679,7 @@ namespace storm { } return rowGroupIndices.get(); } - + template std::vector::index_type> SparseMatrix::swapRowGroupIndices(std::vector&& newRowGrouping) { std::vector result; @@ -689,7 +689,7 @@ namespace storm { } return result; } - + template void SparseMatrix::setRowGroupIndices(std::vector const& newRowGroupIndices) { trivialRowGrouping = false; @@ -700,7 +700,7 @@ namespace storm { bool SparseMatrix::hasTrivialRowGrouping() const { return trivialRowGrouping; } - + template void SparseMatrix::makeRowGroupingTrivial() { if (trivialRowGrouping) { @@ -710,7 +710,7 @@ namespace storm { rowGroupIndices = boost::none; } } - + template storm::storage::BitVector SparseMatrix::getRowFilter(storm::storage::BitVector const& groupConstraint) const { storm::storage::BitVector res(this->getRowCount(), false); @@ -722,7 +722,7 @@ namespace storm { } return res; } - + template storm::storage::BitVector SparseMatrix::getRowFilter(storm::storage::BitVector const& groupConstraint, storm::storage::BitVector const& columnConstraint) const { storm::storage::BitVector result(this->getRowCount(), false); @@ -743,7 +743,7 @@ namespace storm { } return result; } - + template storm::storage::BitVector SparseMatrix::getRowGroupFilter(storm::storage::BitVector const& rowConstraint, bool setIfForAllRowsInGroup) const { STORM_LOG_ASSERT(!this->hasTrivialRowGrouping(), "Tried to get a row group filter but this matrix does not have row groups"); @@ -766,14 +766,14 @@ namespace storm { } return result; } - + template void SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rows) { for (auto row : rows) { makeRowDirac(row, row); } } - + template void SparseMatrix::makeRowGroupsAbsorbing(storm::storage::BitVector const& rowGroupConstraint) { if (!this->hasTrivialRowGrouping()) { @@ -788,19 +788,19 @@ namespace storm { } } } - + template void SparseMatrix::makeRowDirac(index_type row, index_type column) { iterator columnValuePtr = this->begin(row); iterator columnValuePtrEnd = this->end(row); - + // If the row has no elements in it, we cannot make it absorbing, because we would need to move all elements // in the vector of nonzeros otherwise. if (columnValuePtr >= columnValuePtrEnd) { throw storm::exceptions::InvalidStateException() << "Illegal call to SparseMatrix::makeRowDirac: cannot make row " << row << " absorbing, because there is no entry in this row."; } iterator lastColumnValuePtr = this->end(row) - 1; - + // If there is at least one entry in this row, we can set it to one, modify its column value to the // one given by the parameter and set all subsequent elements of this row to zero. // However, we want to preserve that column indices within a row are ascending, so we pick an entry that is close to the desired column index @@ -824,7 +824,7 @@ namespace storm { columnValuePtr->setValue(storm::utility::zero()); } } - + template bool SparseMatrix::compareRows(index_type i1, index_type i2) const { const_iterator end1 = this->end(i1); @@ -841,7 +841,7 @@ namespace storm { } return false; } - + template BitVector SparseMatrix::duplicateRowsInRowgroups() const { BitVector bv(this->getRowCount()); @@ -856,31 +856,31 @@ namespace storm { } return bv; } - + template void SparseMatrix::swapRows(index_type const& row1, index_type const& row2) { if (row1 == row2) { return; } - + // Get the index of the row that has more / less entries than the other. index_type largerRow = getRow(row1).getNumberOfEntries() > getRow(row2).getNumberOfEntries() ? row1 : row2; index_type smallerRow = largerRow == row1 ? row2 : row1; index_type rowSizeDifference = getRow(largerRow).getNumberOfEntries() - getRow(smallerRow).getNumberOfEntries(); - + // Save contents of larger row. auto copyRow = getRow(largerRow); std::vector> largerRowContents(copyRow.begin(), copyRow.end()); - + if (largerRow < smallerRow) { auto writeIt = getRows(largerRow, smallerRow + 1).begin(); - + // Write smaller row to its new position. for (auto& smallerRowEntry : getRow(smallerRow)) { *writeIt = std::move(smallerRowEntry); ++writeIt; } - + // Write the intermediate rows into their correct position. if (!storm::utility::isZero(rowSizeDifference)) { for (auto& intermediateRowEntry : getRows(largerRow + 1, smallerRow)) { @@ -891,15 +891,15 @@ namespace storm { // skip the intermediate rows writeIt = getRow(smallerRow).begin(); } - + // Write the larger row to its new position. for (auto& largerRowEntry : largerRowContents) { *writeIt = std::move(largerRowEntry); ++writeIt; } - + STORM_LOG_ASSERT(writeIt == getRow(smallerRow).end(), "Unexpected position of write iterator."); - + // Update the row indications to account for the shift of indices at where the rows now start. if (!storm::utility::isZero(rowSizeDifference)) { for (index_type row = largerRow + 1; row <= smallerRow; ++row) { @@ -908,14 +908,14 @@ namespace storm { } } else { auto writeIt = getRows(smallerRow, largerRow + 1).end() - 1; - + // Write smaller row to its new position auto copyRow = getRow(smallerRow); for (auto smallerRowEntryIt = copyRow.end() - 1; smallerRowEntryIt != copyRow.begin() - 1; --smallerRowEntryIt) { *writeIt = std::move(*smallerRowEntryIt); --writeIt; } - + // Write the intermediate rows into their correct position. if (!storm::utility::isZero(rowSizeDifference)) { for (auto intermediateRowEntryIt = getRows(smallerRow + 1, largerRow).end() - 1; intermediateRowEntryIt != getRows(smallerRow + 1, largerRow).begin() - 1; --intermediateRowEntryIt) { @@ -926,15 +926,15 @@ namespace storm { // skip the intermediate rows writeIt = getRow(smallerRow).end() - 1; } - + // Write the larger row to its new position. for (auto largerRowEntryIt = largerRowContents.rbegin(); largerRowEntryIt != largerRowContents.rend(); ++largerRowEntryIt) { *writeIt = std::move(*largerRowEntryIt); --writeIt; } - + STORM_LOG_ASSERT(writeIt == getRow(smallerRow).begin() - 1, "Unexpected position of write iterator."); - + // Update row indications. // Update the row indications to account for the shift of indices at where the rows now start. if (!storm::utility::isZero(rowSizeDifference)) { @@ -944,11 +944,11 @@ namespace storm { } } } - + template std::vector SparseMatrix::getRowSumVector() const { std::vector result(this->getRowCount()); - + index_type row = 0; for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++row) { *resultIt = getRowSum(row); @@ -956,7 +956,7 @@ namespace storm { return result; } - + template ValueType SparseMatrix::getConstrainedRowSum(index_type row, storm::storage::BitVector const& constraint) const { ValueType result = storm::utility::zero(); @@ -967,7 +967,7 @@ namespace storm { } return result; } - + template std::vector SparseMatrix::getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const { std::vector result(rowConstraint.getNumberOfSetBits()); @@ -977,7 +977,7 @@ namespace storm { } return result; } - + template std::vector SparseMatrix::getConstrainedRowGroupSumVector(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint) const { std::vector result; @@ -995,7 +995,7 @@ namespace storm { } return result; } - + template SparseMatrix SparseMatrix::getSubmatrix(bool useGroups, storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint, bool insertDiagonalElements, storm::storage::BitVector const& makeZeroColumns) const { if (useGroups) { @@ -1008,14 +1008,14 @@ namespace storm { *it = i; } auto res = getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements, makeZeroColumns); - + // Create a new row grouping that reflects the new sizes of the row groups if the current matrix has a // non trivial row-grouping. if (!this->hasTrivialRowGrouping()) { std::vector newRowGroupIndices; newRowGroupIndices.push_back(0); auto selectedRowIt = rowConstraint.begin(); - + // For this, we need to count how many rows were preserved in every group. for (uint_fast64_t group = 0; group < this->getRowGroupCount(); ++group) { uint_fast64_t newRowCount = 0; @@ -1027,20 +1027,20 @@ namespace storm { newRowGroupIndices.push_back(newRowGroupIndices.back() + newRowCount); } } - + res.trivialRowGrouping = false; res.rowGroupIndices = newRowGroupIndices; } - + return res; } } - + template SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries, storm::storage::BitVector const& makeZeroColumns) const { STORM_LOG_THROW(!rowGroupConstraint.empty() && !columnConstraint.empty(), storm::exceptions::InvalidArgumentException, "Cannot build empty submatrix."); uint_fast64_t submatrixColumnCount = columnConstraint.getNumberOfSetBits(); - + // Start by creating a temporary vector that stores for each index whose bit is set to true the number of // bits that were set before that particular index. std::vector columnBitsSetBeforeIndex = columnConstraint.getNumberOfSetBitsBeforeIndices(); @@ -1049,7 +1049,7 @@ namespace storm { tmp = std::make_unique>(rowGroupConstraint.getNumberOfSetBitsBeforeIndices()); } std::vector const& rowBitsSetBeforeIndex = tmp ? *tmp : columnBitsSetBeforeIndex; - + // Then, we need to determine the number of entries and the number of rows of the submatrix. index_type subEntries = 0; index_type subRows = 0; @@ -1058,17 +1058,17 @@ namespace storm { subRows += rowGroupIndices[index + 1] - rowGroupIndices[index]; for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { bool foundDiagonalElement = false; - + for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) { if (columnConstraint.get(it->getColumn()) && (makeZeroColumns.size() == 0 || !makeZeroColumns.get(it->getColumn())) ) { ++subEntries; - + if (columnBitsSetBeforeIndex[it->getColumn()] == rowBitsSetBeforeIndex[index]) { foundDiagonalElement = true; } } } - + // If requested, we need to reserve one entry more for inserting the diagonal zero entry. if (insertDiagonalEntries && !foundDiagonalElement && rowGroupCount < submatrixColumnCount) { ++subEntries; @@ -1076,10 +1076,10 @@ namespace storm { } ++rowGroupCount; } - + // Create and initialize resulting matrix. SparseMatrixBuilder matrixBuilder(subRows, submatrixColumnCount, subEntries, true, !this->hasTrivialRowGrouping()); - + // Copy over selected entries. rowGroupCount = 0; index_type rowCount = 0; @@ -1090,7 +1090,7 @@ namespace storm { } for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { bool insertedDiagonalElement = false; - + for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) { if (columnConstraint.get(it->getColumn()) && (makeZeroColumns.size() == 0 || !makeZeroColumns.get(it->getColumn()))) { if (columnBitsSetBeforeIndex[it->getColumn()] == rowBitsSetBeforeIndex[index]) { @@ -1110,20 +1110,20 @@ namespace storm { } ++rowGroupCount; } - + return matrixBuilder.build(); } - + template SparseMatrix SparseMatrix::restrictRows(storm::storage::BitVector const& rowsToKeep, bool allowEmptyRowGroups) const { STORM_LOG_ASSERT(rowsToKeep.size() == this->getRowCount(), "Dimensions mismatch."); - + // Count the number of entries of the resulting matrix uint_fast64_t entryCount = 0; for (auto const& row : rowsToKeep) { entryCount += this->getRow(row).getNumberOfEntries(); } - + // Get the smallest row group index such that all row groups with at least this index are empty. uint_fast64_t firstTrailingEmptyRowGroup = this->getRowGroupCount(); for (auto groupIndexIt = this->getRowGroupIndices().rbegin() + 1; groupIndexIt != this->getRowGroupIndices().rend(); ++groupIndexIt) { @@ -1133,7 +1133,7 @@ namespace storm { --firstTrailingEmptyRowGroup; } STORM_LOG_THROW(allowEmptyRowGroups || firstTrailingEmptyRowGroup == this->getRowGroupCount(), storm::exceptions::InvalidArgumentException, "Empty rows are not allowed, but row group " << firstTrailingEmptyRowGroup << " is empty."); - + // build the matrix. The row grouping will always be considered as nontrivial. SparseMatrixBuilder builder(rowsToKeep.getNumberOfSetBits(), this->getColumnCount(), entryCount, true, true, this->getRowGroupCount()); uint_fast64_t newRow = 0; @@ -1150,12 +1150,12 @@ namespace storm { } STORM_LOG_THROW(allowEmptyRowGroups || !rowGroupEmpty, storm::exceptions::InvalidArgumentException, "Empty rows are not allowed, but row group " << rowGroup << " is empty."); } - + // The all remaining row groups will be empty. Note that it is not allowed to call builder.addNewGroup(...) if there are no more rows afterwards. SparseMatrix res = builder.build(); return res; } - + template SparseMatrix SparseMatrix::filterEntries(storm::storage::BitVector const& rowFilter) const { // Count the number of entries in the resulting matrix. @@ -1163,7 +1163,7 @@ namespace storm { for (auto const& row : rowFilter) { entryCount += getRow(row).getNumberOfEntries(); } - + // Build the resulting matrix. SparseMatrixBuilder builder(getRowCount(), getColumnCount(), entryCount); for (auto const& row : rowFilter) { @@ -1172,14 +1172,14 @@ namespace storm { } } SparseMatrix result = builder.build(); - + // Add a row grouping if necessary. if (!hasTrivialRowGrouping()) { result.setRowGroupIndices(getRowGroupIndices()); } return result; } - + template void SparseMatrix::dropZeroEntries() { updateNonzeroEntryCount(); @@ -1200,7 +1200,7 @@ namespace storm { *this = std::move(result); } } - + template SparseMatrix SparseMatrix::selectRowsFromRowGroups(std::vector const& rowGroupToRowIndexMapping, bool insertDiagonalEntries) const { // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for @@ -1209,7 +1209,7 @@ namespace storm { for (index_type rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { // Determine which row we need to select from the current row group. index_type rowToCopy = this->getRowGroupIndices()[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; - + // Iterate through that row and count the number of slots we have to reserve for copying. bool foundDiagonalElement = false; for (const_iterator it = this->begin(rowToCopy), ite = this->end(rowToCopy); it != ite; ++it) { @@ -1222,15 +1222,15 @@ namespace storm { ++subEntries; } } - + // Now create the matrix to be returned with the appropriate size. SparseMatrixBuilder matrixBuilder(rowGroupIndices.get().size() - 1, columnCount, subEntries); - + // Copy over the selected lines from the source matrix. for (index_type rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { // Determine which row we need to select from the current row group. index_type rowToCopy = this->getRowGroupIndices()[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; - + // Iterate through that row and copy the entries. This also inserts a zero element on the diagonal if // there is no entry yet. bool insertedDiagonalElement = false; @@ -1247,7 +1247,7 @@ namespace storm { matrixBuilder.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::zero()); } } - + // Finalize created matrix and return result. return matrixBuilder.build(); } @@ -1317,7 +1317,7 @@ namespace storm { } return result; } - + template SparseMatrix SparseMatrix::transpose(bool joinGroups, bool keepZeros) const { index_type rowCount = this->getColumnCount(); @@ -1329,10 +1329,10 @@ namespace storm { this->updateNonzeroEntryCount(); entryCount = this->getNonzeroEntryCount(); } - + std::vector rowIndications(rowCount + 1); std::vector> columnsAndValues(entryCount); - + // First, we need to count how many entries each column has. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { @@ -1341,17 +1341,17 @@ namespace storm { } } } - + // Now compute the accumulated offsets. for (index_type i = 1; i < rowCount + 1; ++i) { rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; } - + // Create an array that stores the index for the next value to be added for // each row in the transposed matrix. Initially this corresponds to the previously // computed accumulated offsets. std::vector nextIndices = rowIndications; - + // Now we are ready to actually fill in the values of the transposed matrix. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { @@ -1361,17 +1361,17 @@ namespace storm { } } } - + storm::storage::SparseMatrix transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), boost::none); - + return transposedMatrix; } - + template SparseMatrix SparseMatrix::transposeSelectedRowsFromRowGroups(std::vector const& rowGroupChoices, bool keepZeros) const { index_type rowCount = this->getColumnCount(); index_type columnCount = this->getRowGroupCount(); - + // Get the overall entry count as well as the number of entries of each column index_type entryCount = 0; std::vector rowIndications(columnCount + 1); @@ -1384,19 +1384,19 @@ namespace storm { } } } - + // Now compute the accumulated offsets. for (index_type i = 1; i < rowCount + 1; ++i) { rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; } - + std::vector> columnsAndValues(entryCount); - + // Create an array that stores the index for the next value to be added for // each row in the transposed matrix. Initially this corresponds to the previously // computed accumulated offsets. std::vector nextIndices = rowIndications; - + // Now we are ready to actually fill in the values of the transposed matrix. rowGroupChoiceIt = rowGroupChoices.begin(); for (index_type rowGroup = 0; rowGroup < columnCount; ++rowGroup, ++rowGroupChoiceIt) { @@ -1407,16 +1407,16 @@ namespace storm { } } } - + return storm::storage::SparseMatrix(std::move(columnCount), std::move(rowIndications), std::move(columnsAndValues), boost::none); } - + template void SparseMatrix::convertToEquationSystem() { invertDiagonal(); negateAllNonDiagonalEntries(); } - + template void SparseMatrix::invertDiagonal() { // Now iterate over all row groups and set the diagonal elements to the inverted value. @@ -1439,14 +1439,14 @@ namespace storm { foundDiagonalElement = true; } } - + // Throw an exception if a row did not have an element on the diagonal. if (!foundDiagonalElement) { throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::invertDiagonal: matrix is missing diagonal entries."; } } } - + template void SparseMatrix::negateAllNonDiagonalEntries() { // Iterate over all row groups and negate all the elements that are not on the diagonal. @@ -1458,7 +1458,7 @@ namespace storm { } } } - + template void SparseMatrix::deleteDiagonalEntries() { // Iterate over all rows and negate all the elements that are not on the diagonal. @@ -1471,15 +1471,15 @@ namespace storm { } } } - + template typename std::pair, std::vector> SparseMatrix::getJacobiDecomposition() const { STORM_LOG_THROW(this->getRowCount() == this->getColumnCount(), storm::exceptions::InvalidArgumentException, "Canno compute Jacobi decomposition of non-square matrix."); - + // Prepare the resulting data structures. SparseMatrixBuilder luBuilder(this->getRowCount(), this->getColumnCount()); std::vector invertedDiagonal(rowCount); - + // Copy entries to the appropriate matrices. for (index_type rowNumber = 0; rowNumber < rowCount; ++rowNumber) { for (const_iterator it = this->begin(rowNumber), ite = this->end(rowNumber); it != ite; ++it) { @@ -1490,22 +1490,22 @@ namespace storm { } } } - + return std::make_pair(luBuilder.build(), std::move(invertedDiagonal)); } - + #ifdef STORM_HAVE_CARL template<> typename std::pair, std::vector> SparseMatrix::getJacobiDecomposition() const { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); } - + template<> typename std::pair, std::vector> SparseMatrix::getJacobiDecomposition() const { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); } #endif - + template template ResultValueType SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, index_type const& row) const { @@ -1513,7 +1513,7 @@ namespace storm { typename storm::storage::SparseMatrix::const_iterator ite1 = this->end(row); typename storm::storage::SparseMatrix::const_iterator it2 = otherMatrix.begin(row); typename storm::storage::SparseMatrix::const_iterator ite2 = otherMatrix.end(row); - + ResultValueType result = storm::utility::zero(); for (;it1 != ite1 && it2 != ite2; ++it1) { if (it1->getColumn() < it2->getColumn()) { @@ -1529,7 +1529,7 @@ namespace storm { } return result; } - + template template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const { @@ -1540,7 +1540,7 @@ namespace storm { } return result; } - + template void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand) const { // If the vector and the result are aliases and this is not set to be allowed, we need and temporary vector. @@ -1553,14 +1553,14 @@ namespace storm { } else { target = &result; } - + this->multiplyWithVectorForward(vector, *target, summand); - + if (target == &temporary) { std::swap(result, *target); } } - + template void SparseMatrix::multiplyWithVectorForward(std::vector const& vector, std::vector& result, std::vector const* summand) const { const_iterator it = this->begin(); @@ -1572,7 +1572,7 @@ namespace storm { if (summand) { summandIterator = summand->begin(); } - + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++summandIterator) { ValueType newValue; if (summand) { @@ -1580,15 +1580,15 @@ namespace storm { } else { newValue = storm::utility::zero(); } - + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { newValue += it->getValue() * vector[it->getColumn()]; } - + *resultIterator = newValue; } } - + template void SparseMatrix::multiplyWithVectorBackward(std::vector const& vector, std::vector& result, std::vector const* summand) const { const_iterator it = this->end() - 1; @@ -1600,7 +1600,7 @@ namespace storm { if (summand) { summandIterator = summand->end() - 1; } - + for (; resultIterator != resultIteratorEnd; --rowIterator, --resultIterator, --summandIterator) { ValueType newValue; if (summand) { @@ -1608,15 +1608,15 @@ namespace storm { } else { newValue = storm::utility::zero(); } - + for (ite = this->begin() + *rowIterator - 1; it != ite; --it) { newValue += (it->getValue() * vector[it->getColumn()]); } - + *resultIterator = newValue; } } - + #ifdef STORM_HAVE_INTELTBB template class TbbMultAddFunctor { @@ -1624,11 +1624,11 @@ namespace storm { typedef typename storm::storage::SparseMatrix::index_type index_type; typedef typename storm::storage::SparseMatrix::value_type value_type; typedef typename storm::storage::SparseMatrix::const_iterator const_iterator; - + TbbMultAddFunctor(std::vector> const& columnsAndEntries, std::vector const& rowIndications, std::vector const& x, std::vector& result, std::vector const* summand) : columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand) { // Intentionally left empty. } - + void operator()(tbb::blocked_range const& range) const { index_type startRow = range.begin(); index_type endRow = range.end(); @@ -1641,18 +1641,18 @@ namespace storm { if (summand) { summandIterator = summand->begin() + startRow; } - + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++summandIterator) { ValueType newValue = summand ? *summandIterator : storm::utility::zero(); - + for (ite = columnsAndEntries.begin() + *(rowIterator + 1); it != ite; ++it) { newValue += it->getValue() * x[it->getColumn()]; } - + *resultIterator = newValue; } } - + private: std::vector> const& columnsAndEntries; std::vector const& rowIndications; @@ -1660,7 +1660,7 @@ namespace storm { std::vector& result; std::vector const* summand; }; - + template void SparseMatrix::multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand) const { if (&vector == &result) { @@ -1673,7 +1673,7 @@ namespace storm { } } #endif - + template ValueType SparseMatrix::multiplyRowWithVector(index_type row, std::vector const& vector) const { ValueType result = storm::utility::zero(); @@ -1683,7 +1683,7 @@ namespace storm { } return result; } - + template void SparseMatrix::performSuccessiveOverRelaxationStep(ValueType omega, std::vector& x, std::vector const& b) const { const_iterator it = this->end() - 1; @@ -1692,13 +1692,13 @@ namespace storm { typename std::vector::const_iterator bIt = b.end() - 1; typename std::vector::iterator resultIterator = x.end() - 1; typename std::vector::iterator resultIteratorEnd = x.begin() - 1; - + index_type currentRow = getRowCount(); for (; resultIterator != resultIteratorEnd; --rowIterator, --resultIterator, --bIt) { --currentRow; ValueType tmpValue = storm::utility::zero(); ValueType diagonalElement = storm::utility::zero(); - + for (ite = this->begin() + *rowIterator - 1; it != ite; --it) { if (it->getColumn() != currentRow) { tmpValue += it->getValue() * x[it->getColumn()]; @@ -1710,14 +1710,14 @@ namespace storm { *resultIterator = ((storm::utility::one() - omega) * *resultIterator) + (omega / diagonalElement) * (*bIt - tmpValue); } } - + #ifdef STORM_HAVE_CARL template<> void SparseMatrix::performSuccessiveOverRelaxationStep(Interval, std::vector&, std::vector const&) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif - + template void SparseMatrix::performWalkerChaeStep(std::vector const& x, std::vector const& columnSums, std::vector const& b, std::vector const& ax, std::vector& result) const { const_iterator it = this->begin(); @@ -1735,7 +1735,7 @@ namespace storm { result[it->getColumn()] += it->getValue() * (b[row] / ax[row]); } } - + auto xIterator = x.begin(); auto sumsIterator = columnSums.begin(); for (auto& entry : result) { @@ -1751,7 +1751,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif - + template void SparseMatrix::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { if (dir == OptimizationDirection::Minimize) { @@ -1776,36 +1776,36 @@ namespace storm { if (choices) { choiceIt = choices->begin(); } - + // Variables for correctly tracking choices (only update if new choice is strictly better). ValueType oldSelectedChoiceValue; uint64_t selectedChoice; - + uint64_t currentRow = 0; for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) { ValueType currentValue = storm::utility::zero(); - + // Only multiply and reduce if there is at least one row in the group. if (*rowGroupIt < *(rowGroupIt + 1)) { if (summand) { currentValue = *summandIt; ++summandIt; } - + for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { currentValue += elementIt->getValue() * vector[elementIt->getColumn()]; } - + if (choices) { selectedChoice = 0; if (*choiceIt == 0) { oldSelectedChoiceValue = currentValue; } } - + ++rowIt; ++currentRow; - + for (; currentRow < *(rowGroupIt + 1); ++rowIt, ++currentRow) { ValueType newValue = summand ? *summandIt : storm::utility::zero(); for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { @@ -1815,7 +1815,7 @@ namespace storm { if (choices && currentRow == *choiceIt + *rowGroupIt) { oldSelectedChoiceValue = newValue; } - + if (compare(newValue, currentValue)) { currentValue = newValue; if (choices) { @@ -1826,7 +1826,7 @@ namespace storm { ++summandIt; } } - + // Finally write value to target vector. *resultIt = currentValue; if (choices && compare(currentValue, oldSelectedChoiceValue)) { @@ -1835,14 +1835,14 @@ namespace storm { } } } - + #ifdef STORM_HAVE_CARL template<> void SparseMatrix::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif - + template void SparseMatrix::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { if (dir == storm::OptimizationDirection::Minimize) { @@ -1851,7 +1851,7 @@ namespace storm { multiplyAndReduceBackward>(rowGroupIndices, vector, summand, result, choices); } } - + template template void SparseMatrix::multiplyAndReduceBackward(std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { @@ -1867,22 +1867,22 @@ namespace storm { if (choices) { choiceIt = choices->end() - 1; } - + // Variables for correctly tracking choices (only update if new choice is strictly better). ValueType oldSelectedChoiceValue; uint64_t selectedChoice; - + uint64_t currentRow = this->getRowCount() - 1; for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt) { ValueType currentValue = storm::utility::zero(); - + // Only multiply and reduce if there is at least one row in the group. if (*rowGroupIt < *(rowGroupIt + 1)) { if (summand) { currentValue = *summandIt; --summandIt; } - + for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) { currentValue += elementIt->getValue() * vector[elementIt->getColumn()]; } @@ -1894,17 +1894,17 @@ namespace storm { } --rowIt; --currentRow; - + for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, --currentRow, ++i, --summandIt) { ValueType newValue = summand ? *summandIt : storm::utility::zero(); for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) { newValue += elementIt->getValue() * vector[elementIt->getColumn()]; } - + if (choices && currentRow == *choiceIt + *rowGroupIt) { oldSelectedChoiceValue = newValue; } - + if (compare(newValue, currentValue)) { currentValue = newValue; if (choices) { @@ -1912,7 +1912,7 @@ namespace storm { } } } - + // Finally write value to target vector. *resultIt = currentValue; if (choices && compare(currentValue, oldSelectedChoiceValue)) { @@ -1921,14 +1921,14 @@ namespace storm { } } } - + #ifdef STORM_HAVE_CARL template<> void SparseMatrix::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif - + #ifdef STORM_HAVE_INTELTBB template class TbbMultAddReduceFunctor { @@ -1936,16 +1936,16 @@ namespace storm { typedef typename storm::storage::SparseMatrix::index_type index_type; typedef typename storm::storage::SparseMatrix::value_type value_type; typedef typename storm::storage::SparseMatrix::const_iterator const_iterator; - + TbbMultAddReduceFunctor(std::vector const& rowGroupIndices, std::vector> const& columnsAndEntries, std::vector const& rowIndications, std::vector const& x, std::vector& result, std::vector const* summand, std::vector* choices) : rowGroupIndices(rowGroupIndices), columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand), choices(choices) { // Intentionally left empty. } - + void operator()(tbb::blocked_range const& range) const { auto groupIt = rowGroupIndices.begin() + range.begin(); auto groupIte = rowGroupIndices.begin() + range.end(); - + auto rowIt = rowIndications.begin() + *groupIt; auto elementIt = columnsAndEntries.begin() + *rowIt; typename std::vector::const_iterator summandIt; @@ -1956,9 +1956,9 @@ namespace storm { if (choices) { choiceIt = choices->begin() + range.begin(); } - + auto resultIt = result.begin() + range.begin(); - + // Variables for correctly tracking choices (only update if new choice is strictly better). ValueType oldSelectedChoiceValue; uint64_t selectedChoice; @@ -1966,7 +1966,7 @@ namespace storm { uint64_t currentRow = *groupIt; for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { ValueType currentValue = storm::utility::zero(); - + // Only multiply and reduce if there is at least one row in the group. if (*groupIt < *(groupIt + 1)) { if (summand) { @@ -1977,23 +1977,23 @@ namespace storm { for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { currentValue += elementIt->getValue() * x[elementIt->getColumn()]; } - + if (choices) { selectedChoice = 0; if (*choiceIt == 0) { oldSelectedChoiceValue = currentValue; } } - + ++rowIt; ++currentRow; - + for (; currentRow < *(groupIt + 1); ++rowIt, ++currentRow, ++summandIt) { ValueType newValue = summand ? *summandIt : storm::utility::zero(); for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { newValue += elementIt->getValue() * x[elementIt->getColumn()]; } - + if (choices && currentRow == *choiceIt + *groupIt) { oldSelectedChoiceValue = newValue; } @@ -2005,7 +2005,7 @@ namespace storm { } } } - + // Finally write value to target vector. *resultIt = currentValue; if (choices && compare(currentValue, oldSelectedChoiceValue)) { @@ -2014,7 +2014,7 @@ namespace storm { } } } - + private: Compare compare; std::vector const& rowGroupIndices; @@ -2025,7 +2025,7 @@ namespace storm { std::vector const* summand; std::vector* choices; }; - + template void SparseMatrix::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { if (dir == storm::OptimizationDirection::Minimize) { @@ -2034,7 +2034,7 @@ namespace storm { tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices)); } } - + #ifdef STORM_HAVE_CARL template<> void SparseMatrix::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { @@ -2042,10 +2042,10 @@ namespace storm { } #endif #endif - + template void SparseMatrix::multiplyAndReduce(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { - + // If the vector and the result are aliases, we need and temporary vector. std::vector* target; std::vector temporary; @@ -2056,21 +2056,21 @@ namespace storm { } else { target = &result; } - + this->multiplyAndReduceForward(dir, rowGroupIndices, vector, summand, *target, choices); if (target == &temporary) { std::swap(temporary, result); } } - + template void SparseMatrix::multiplyVectorWithMatrix(std::vector const& vector, std::vector& result) const { const_iterator it = this->begin(); const_iterator ite; std::vector::const_iterator rowIterator = rowIndications.begin(); std::vector::const_iterator rowIteratorEnd = rowIndications.end(); - + uint_fast64_t currentRow = 0; for (; rowIterator != rowIteratorEnd - 1; ++rowIterator) { for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { @@ -2079,7 +2079,7 @@ namespace storm { ++currentRow; } } - + template void SparseMatrix::scaleRowsInPlace(std::vector const& factors) { STORM_LOG_ASSERT(factors.size() == this->getRowCount(), "Can not scale rows: Number of rows and number of scaling factors do not match."); @@ -2091,7 +2091,7 @@ namespace storm { ++row; } } - + template void SparseMatrix::divideRowsInPlace(std::vector const& divisors) { STORM_LOG_ASSERT(divisors.size() == this->getRowCount(), "Can not divide rows: Number of rows and number of divisors do not match."); @@ -2104,34 +2104,34 @@ namespace storm { ++row; } } - + #ifdef STORM_HAVE_CARL template<> void SparseMatrix::divideRowsInPlace(std::vector const&) { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); } #endif - + template typename SparseMatrix::const_rows SparseMatrix::getRows(index_type startRow, index_type endRow) const { return const_rows(this->columnsAndValues.begin() + this->rowIndications[startRow], this->rowIndications[endRow] - this->rowIndications[startRow]); } - + template typename SparseMatrix::rows SparseMatrix::getRows(index_type startRow, index_type endRow) { return rows(this->columnsAndValues.begin() + this->rowIndications[startRow], this->rowIndications[endRow] - this->rowIndications[startRow]); } - + template typename SparseMatrix::const_rows SparseMatrix::getRow(index_type row) const { return getRows(row, row + 1); } - + template typename SparseMatrix::rows SparseMatrix::getRow(index_type row) { return getRows(row, row + 1); } - + template typename SparseMatrix::const_rows SparseMatrix::getRow(index_type rowGroup, index_type offset) const { STORM_LOG_ASSERT(rowGroup < this->getRowGroupCount(), "Row group is out-of-bounds."); @@ -2142,7 +2142,7 @@ namespace storm { return getRow(this->getRowGroupIndices()[rowGroup] + offset); } } - + template typename SparseMatrix::rows SparseMatrix::getRow(index_type rowGroup, index_type offset) { STORM_LOG_ASSERT(rowGroup < this->getRowGroupCount(), "Row group is out-of-bounds."); @@ -2154,8 +2154,8 @@ namespace storm { return getRow(rowGroup + offset); } } - - + + template typename SparseMatrix::const_rows SparseMatrix::getRowGroup(index_type rowGroup) const { STORM_LOG_ASSERT(rowGroup < this->getRowGroupCount(), "Row group is out-of-bounds."); @@ -2165,7 +2165,7 @@ namespace storm { return getRows(rowGroup, rowGroup + 1); } } - + template typename SparseMatrix::rows SparseMatrix::getRowGroup(index_type rowGroup) { STORM_LOG_ASSERT(rowGroup < this->getRowGroupCount(), "Row group is out-of-bounds."); @@ -2175,32 +2175,32 @@ namespace storm { return getRows(rowGroup, rowGroup + 1); } } - + template typename SparseMatrix::const_iterator SparseMatrix::begin(index_type row) const { return this->columnsAndValues.begin() + this->rowIndications[row]; } - + template typename SparseMatrix::iterator SparseMatrix::begin(index_type row) { return this->columnsAndValues.begin() + this->rowIndications[row]; } - + template typename SparseMatrix::const_iterator SparseMatrix::end(index_type row) const { return this->columnsAndValues.begin() + this->rowIndications[row + 1]; } - + template typename SparseMatrix::iterator SparseMatrix::end(index_type row) { return this->columnsAndValues.begin() + this->rowIndications[row + 1]; } - + template typename SparseMatrix::const_iterator SparseMatrix::end() const { return this->columnsAndValues.begin() + this->rowIndications[rowCount]; } - + template typename SparseMatrix::iterator SparseMatrix::end() { return this->columnsAndValues.begin() + this->rowIndications[rowCount]; @@ -2214,7 +2214,7 @@ namespace storm { } return sum; } - + template typename SparseMatrix::index_type SparseMatrix::getNonconstantEntryCount() const { index_type nonConstEntries = 0; @@ -2225,7 +2225,7 @@ namespace storm { } return nonConstEntries; } - + template typename SparseMatrix::index_type SparseMatrix::getNonconstantRowGroupCount() const { index_type nonConstRowGroups = 0; @@ -2239,7 +2239,7 @@ namespace storm { } return nonConstRowGroups; } - + template bool SparseMatrix::isProbabilistic() const { storm::utility::ConstantsComparator comparator; @@ -2258,7 +2258,7 @@ namespace storm { } return true; } - + template template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const { @@ -2269,7 +2269,7 @@ namespace storm { if (!this->hasTrivialRowGrouping() && matrix.hasTrivialRowGrouping()) return false; if (!this->hasTrivialRowGrouping() && !matrix.hasTrivialRowGrouping() && this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false; if (this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false; - + // Check the subset property for all rows individually. for (index_type row = 0; row < this->getRowCount(); ++row) { auto it2 = matrix.begin(row); @@ -2286,7 +2286,7 @@ namespace storm { } return true; } - + template bool SparseMatrix::isIdentityMatrix() const { if (this->getRowCount() != this->getColumnCount()) { @@ -2326,7 +2326,7 @@ namespace storm { return result; } - + template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix) { // Print column numbers in header. @@ -2335,22 +2335,22 @@ namespace storm { out << i << "\t"; } out << std::endl; - + // Iterate over all row groups. for (typename SparseMatrix::index_type group = 0; group < matrix.getRowGroupCount(); ++group) { out << "\t---- group " << group << "/" << (matrix.getRowGroupCount() - 1) << " ---- " << std::endl; typename SparseMatrix::index_type start = matrix.hasTrivialRowGrouping() ? group : matrix.getRowGroupIndices()[group]; typename SparseMatrix::index_type end = matrix.hasTrivialRowGrouping() ? group + 1 : matrix.getRowGroupIndices()[group + 1]; - + for (typename SparseMatrix::index_type i = start; i < end; ++i) { typename SparseMatrix::index_type nextIndex = matrix.rowIndications[i]; - + // Print the actual row. out << i << "\t(\t"; typename SparseMatrix::index_type currentRealIndex = 0; while (currentRealIndex < matrix.columnCount) { if (nextIndex < matrix.rowIndications[i + 1] && currentRealIndex == matrix.columnsAndValues[nextIndex].getColumn()) { - out << matrix.columnsAndValues[nextIndex].getValue() << "\t"; + out << std::setprecision(3) << matrix.columnsAndValues[nextIndex].getValue() << "\t"; ++nextIndex; } else { out << "0\t"; @@ -2360,17 +2360,17 @@ namespace storm { out << "\t)\t" << i << std::endl; } } - + // Print column numbers in footer. out << "\t\t"; for (typename SparseMatrix::index_type i = 0; i < matrix.getColumnCount(); ++i) { out << i << "\t"; } out << std::endl; - + return out; } - + template void SparseMatrix::printAsMatlabMatrix(std::ostream& out) const { // Iterate over all row groups. @@ -2378,7 +2378,7 @@ namespace storm { STORM_LOG_ASSERT(this->getRowGroupSize(group) == 1, "Incorrect row group size."); for (typename SparseMatrix::index_type i = this->getRowGroupIndices()[group]; i < this->getRowGroupIndices()[group + 1]; ++i) { typename SparseMatrix::index_type nextIndex = this->rowIndications[i]; - + // Print the actual row. out << i << "\t("; typename SparseMatrix::index_type currentRealIndex = 0; @@ -2395,11 +2395,11 @@ namespace storm { } } } - + template std::size_t SparseMatrix::hash() const { std::size_t result = 0; - + boost::hash_combine(result, this->getRowCount()); boost::hash_combine(result, this->getColumnCount()); boost::hash_combine(result, this->getEntryCount()); @@ -2408,11 +2408,11 @@ namespace storm { if (!this->hasTrivialRowGrouping()) { boost::hash_combine(result, boost::hash_range(rowGroupIndices.get().begin(), rowGroupIndices.get().end())); } - + return result; } - - + + #ifdef STORM_HAVE_CARL std::set getVariables(SparseMatrix const& matrix) { @@ -2422,9 +2422,9 @@ namespace storm { } return result; } - + #endif - + // Explicitly instantiate the entry, builder and the matrix. // double template class MatrixEntry::index_type, double>; @@ -2435,10 +2435,10 @@ namespace storm { template double SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + template class MatrixEntry; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); - + // float template class MatrixEntry::index_type, float>; template std::ostream& operator<<(std::ostream& out, MatrixEntry::index_type, float> const& entry); @@ -2448,7 +2448,7 @@ namespace storm { template float SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + // int template class MatrixEntry::index_type, int>; template std::ostream& operator<<(std::ostream& out, MatrixEntry::index_type, int> const& entry); @@ -2456,7 +2456,7 @@ namespace storm { template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + // state_type template class MatrixEntry::index_type, storm::storage::sparse::state_type>; template std::ostream& operator<<(std::ostream& out, MatrixEntry::index_type, storm::storage::sparse::state_type> const& entry); @@ -2464,10 +2464,10 @@ namespace storm { template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + #ifdef STORM_HAVE_CARL // Rational Numbers - + #if defined(STORM_HAVE_CLN) template class MatrixEntry::index_type, ClnRationalNumber>; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -2478,7 +2478,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif - + #if defined(STORM_HAVE_GMP) template class MatrixEntry::index_type, GmpRationalNumber>; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -2489,7 +2489,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif - + // Rational Function template class MatrixEntry::index_type, RationalFunction>; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -2505,7 +2505,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + // Intervals template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template class MatrixEntry::index_type, Interval>; @@ -2515,13 +2515,10 @@ namespace storm { template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif - - - } // namespace storage -} // namespace storm - + } // namespace storage +} // namespace storm