proper EC elimination in hybrid helper

src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp

@@ -19,6 +19,7 @@ namespace storm {
 
             template<typename ValueType>
             ValueType BaierUpperRewardBoundsComputer<ValueType>::computeUpperBound() {
+                STORM_LOG_TRACE("Computing upper reward bounds using variant-2 of Baier et al.");
                 std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
                 
                 std::vector<uint64_t> stateToScc(transitionMatrix.getRowGroupCount());

src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp

@@ -23,6 +23,7 @@ namespace storm {
             
             template<typename ValueType>
             std::vector<ValueType> DsMpiDtmcUpperRewardBoundsComputer<ValueType>::computeUpperBounds() {
+                STORM_LOG_TRACE("Computing upper reward bounds using DS-MPI.");
                 std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();
                 sweep();
                 ValueType lambda = computeLambda();

src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp

@@ -14,6 +14,8 @@
 #include "storm/models/symbolic/StandardRewardModel.h"
 
 #include "storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h"
+#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h"
+#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h"
 #include "storm/modelchecker/results/SymbolicQualitativeCheckResult.h"
 #include "storm/modelchecker/results/SymbolicQuantitativeCheckResult.h"
 #include "storm/modelchecker/results/HybridQuantitativeCheckResult.h"
@@ -32,6 +34,7 @@ namespace storm {
                 boost::optional<SparseMdpEndComponentInformation<ValueType>> ecInformation;
                 boost::optional<std::vector<uint64_t>> initialScheduler;
                 storm::storage::BitVector properMaybeStates;
+                boost::optional<std::vector<ValueType>> oneStepTargetProbabilities;
             };
             
             template <typename ValueType>
@@ -66,15 +69,17 @@ namespace storm {
             }
             
             template <typename ValueType>
-            void eliminateExtendedStatesFromExplicitRepresentation(std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>>& explicitRepresentation, std::vector<uint64_t>& scheduler, storm::storage::BitVector const& properMaybeStates) {
+            void eliminateExtendedStatesFromExplicitRepresentation(std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>>& explicitRepresentation, boost::optional<std::vector<uint64_t>>& scheduler, storm::storage::BitVector const& properMaybeStates) {
+                if (scheduler) {
                     // Eliminate superfluous entries from the scheduler.
                     uint64_t position = 0;
                     for (auto state : properMaybeStates) {
-                    scheduler[position] = scheduler[state];
+                        scheduler.get()[position] = scheduler.get()[state];
                         position++;
                     }
-                scheduler.resize(properMaybeStates.getNumberOfSetBits());
-                scheduler.shrink_to_fit();
+                    scheduler.get().resize(properMaybeStates.getNumberOfSetBits());
+                    scheduler.get().shrink_to_fit();
+                }
                 
                 // Treat the matrix.
                 explicitRepresentation.first = explicitRepresentation.first.getSubmatrix(true, properMaybeStates, properMaybeStates);
@@ -103,7 +108,7 @@ namespace storm {
                     oneStepProbabilities = std::move(subvector);
                 } else {
                     STORM_LOG_DEBUG("Not eliminating ECs as there are none.");
-                    eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates);
+                    eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates);
                     oneStepProbabilities = explicitRepresentation.first.getConstrainedRowGroupSumVector(solverRequirementsData.properMaybeStates, targetStates);
                 }
             }
@@ -192,14 +197,11 @@ namespace storm {
                             storm::dd::Add<DdType, ValueType> subvector = submatrix * prob1StatesAsColumn;
                             subvector = subvector.sumAbstract(model.getColumnVariables());
                             
-                            // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                            std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd<uint_fast64_t>().sumAbstract(model.getNondeterminismVariables()).toVector(odd);
-                            
                             // Finally cut away all columns targeting non-maybe states.
                             submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
                             
                             // Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
-                            explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                            explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd);
 
                             if (requirements.requiresValidInitialScheduler()) {
                                 solverRequirementsData.initialScheduler = computeValidInitialSchedulerForUntilProbabilities<ValueType>(explicitRepresentation.first, explicitRepresentation.second);
@@ -280,9 +282,6 @@ namespace storm {
                     storm::dd::Add<DdType, ValueType> prob1StatesAsColumn = psiStates.template toAdd<ValueType>().swapVariables(model.getRowColumnMetaVariablePairs());
                     storm::dd::Add<DdType, ValueType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
                     
-                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd<uint_fast64_t>().sumAbstract(model.getNondeterminismVariables()).toVector(odd);
-                    
                     // Finally cut away all columns targeting non-maybe states.
                     submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
                     
@@ -290,7 +289,7 @@ namespace storm {
                     std::vector<ValueType> x(maybeStates.getNonZeroCount(), storm::utility::zero<ValueType>());
                     
                     // Translate the symbolic matrix/vector to their explicit representations.
-                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd);
                     
                     std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
                     solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound);
@@ -338,12 +337,8 @@ namespace storm {
                 // Create the solution vector.
                 std::vector<ValueType> x(model.getNumberOfStates(), storm::utility::zero<ValueType>());
                 
-                // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                storm::dd::Add<DdType, uint_fast64_t> stateActionAdd = (transitionMatrix.notZero().existsAbstract(model.getColumnVariables()) || totalRewardVector.notZero()).template toAdd<uint_fast64_t>();
-                std::vector<uint_fast64_t> rowGroupSizes = stateActionAdd.sumAbstract(model.getNondeterminismVariables()).toVector(odd);
-                
                 // Translate the symbolic matrix/vector to their explicit representations.
-                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = transitionMatrix.toMatrixVector(totalRewardVector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = transitionMatrix.toMatrixVector(totalRewardVector, model.getNondeterminismVariables(), odd, odd);
                 
                 // Perform the matrix-vector multiplication.
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
@@ -386,7 +381,12 @@ namespace storm {
             }
             
             template<typename ValueType>
-            void eliminateEndComponentsAndTargetStatesReachabilityRewards(std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>>& explicitRepresentation, SolverRequirementsData<ValueType>& solverRequirementsData) {
+            std::vector<ValueType> computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(storm::storage::SparseMatrix<ValueType> const& extendedMatrix, storm::storage::BitVector const& properMaybeStates, storm::storage::BitVector const& targetStates) {
+                return extendedMatrix.getConstrainedRowGroupSumVector(properMaybeStates, targetStates);
+            }
+            
+            template<typename ValueType>
+            void eliminateEndComponentsAndTargetStatesReachabilityRewards(std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>>& explicitRepresentation, SolverRequirementsData<ValueType>& solverRequirementsData, storm::storage::BitVector const& targetStates, bool computeOneStepTargetProbabilities) {
 
                 // Get easier handles to the data.
                 auto& transitionMatrix = explicitRepresentation.first;
@@ -431,12 +431,31 @@ namespace storm {
                 // Only do more work if there are actually end-components.
                 if (doDecomposition && !endComponentDecomposition.empty()) {
                     STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " EC(s).");
-                    std::vector<ValueType> subvector;
-                    SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, rewardVector, solverRequirementsData.properMaybeStates, transitionMatrix, subvector);
+                    if (computeOneStepTargetProbabilities) {
+                        solverRequirementsData.oneStepTargetProbabilities = std::vector<ValueType>(solverRequirementsData.properMaybeStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
+                    }
+                    
+                    std::vector<ValueType> subvector(solverRequirementsData.properMaybeStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
+                    SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, solverRequirementsData.properMaybeStates, computeOneStepTargetProbabilities ? &targetStates : nullptr, nullptr, &rewardVector, transitionMatrix, computeOneStepTargetProbabilities ? &solverRequirementsData.oneStepTargetProbabilities.get() : nullptr, &subvector);
                     rewardVector = std::move(subvector);
                 } else {
                     STORM_LOG_DEBUG("Not eliminating ECs as there are none.");
-                    eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates);
+                    if (computeOneStepTargetProbabilities) {
+                        solverRequirementsData.oneStepTargetProbabilities = computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates);
+                    }
+                    eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates);
+                }
+            }
+            
+            template<typename ValueType>
+            void setUpperRewardBounds(storm::solver::MinMaxLinearEquationSolver<ValueType>& solver, storm::OptimizationDirection const& direction, storm::storage::SparseMatrix<ValueType> const& submatrix, std::vector<ValueType> const& choiceRewards, std::vector<ValueType> const& oneStepTargetProbabilities) {
+                // For the min-case, we use DS-MPI, for the max-case variant 2 of the Baier et al. paper (CAV'17).
+                if (direction == storm::OptimizationDirection::Minimize) {
+                    DsMpiMdpUpperRewardBoundsComputer<ValueType> dsmpi(submatrix, choiceRewards, oneStepTargetProbabilities);
+                    solver.setUpperBounds(dsmpi.computeUpperBounds());
+                } else {
+                    BaierUpperRewardBoundsComputer<ValueType> baier(submatrix, choiceRewards, oneStepTargetProbabilities);
+                    solver.setUpperBound(baier.computeUpperBound());
                 }
             }
             
@@ -484,6 +503,11 @@ namespace storm {
                                 clearedRequirements.clearValidInitialScheduler();
                             }
                             clearedRequirements.clearLowerBounds();
+                            if (clearedRequirements.requiresUpperBounds()) {
+                                STORM_LOG_DEBUG("Computing upper bounds, because the solver requires it.");
+                                extendMaybeStates = true;
+                                clearedRequirements.clearUpperBounds();
+                            }
                             STORM_LOG_THROW(clearedRequirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
                         }
 
@@ -500,42 +524,39 @@ namespace storm {
                         // (a) transitions from non-maybe states, and
                         // (b) the choices in the transition matrix that lead to a state that is neither a maybe state
                         // nor a target state ('infinity choices').
-                        storm::dd::Add<DdType, ValueType> choiceFilterAdd = (transitionMatrixBdd && maybeStatesWithTargetStates.renameVariables(model.getRowVariables(), model.getColumnVariables())).existsAbstract(model.getColumnVariables()).template toAdd<ValueType>();
-                        storm::dd::Add<DdType, ValueType> submatrix = transitionMatrix * maybeStatesAdd * choiceFilterAdd;
+                        storm::dd::Add<DdType, ValueType> choiceFilterAdd = maybeStatesAdd * (transitionMatrixBdd && maybeStatesWithTargetStates.renameVariables(model.getRowVariables(), model.getColumnVariables())).existsAbstract(model.getColumnVariables()).template toAdd<ValueType>();
+                        storm::dd::Add<DdType, ValueType> submatrix = transitionMatrix * choiceFilterAdd;
                         
                         // Then compute the reward vector to use in the computation.
-                        storm::dd::Add<DdType, ValueType> subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables());
-                        if (!rewardModel.hasStateActionRewards() && !rewardModel.hasTransitionRewards()) {
-                            // If the reward model neither has state-action nor transition rewards, we need to multiply
-                            // it with the legal nondetermism encodings in each state.
-                            subvector *= choiceFilterAdd;
-                        }
-                        
-                        // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                        storm::dd::Add<DdType, uint_fast64_t> stateActionAdd = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd<uint_fast64_t>();
-                        std::vector<uint_fast64_t> rowGroupSizes = stateActionAdd.sumAbstract(model.getNondeterminismVariables()).toVector(odd);
+                        storm::dd::Add<DdType, ValueType> subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, choiceFilterAdd, submatrix, model.getColumnVariables());
 
                         // Finally cut away all columns targeting non-maybe states (or non-(maybe or target) states, respectively).
                         submatrix *= extendMaybeStates ? maybeStatesWithTargetStates.swapVariables(model.getRowColumnMetaVariablePairs()).template toAdd<ValueType>() : maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
                         
                         // Translate the symbolic matrix/vector to their explicit representations.
-                        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd);
                         
                         // Fulfill the solver's requirements.
                         SolverRequirementsData<ValueType> solverRequirementsData;
-                        if (requirements.requiresNoEndComponents() || requirements.requiresValidInitialScheduler()) {
+                        if (extendMaybeStates) {
                             storm::storage::BitVector targetStates = computeTargetStatesForReachabilityRewardsFromExplicitRepresentation(explicitRepresentation.first);
                             solverRequirementsData.properMaybeStates = ~targetStates;
 
                             if (requirements.requiresNoEndComponents()) {
-                                eliminateEndComponentsAndTargetStatesReachabilityRewards(explicitRepresentation, solverRequirementsData);
+                                eliminateEndComponentsAndTargetStatesReachabilityRewards(explicitRepresentation, solverRequirementsData, targetStates, requirements.requiresUpperBounds());
                             } else {
+                                if (requirements.requiresValidInitialScheduler()) {
                                     // Compute a valid initial scheduler.
                                     solverRequirementsData.initialScheduler = computeValidInitialSchedulerForReachabilityRewards<ValueType>(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates);
+                                }
+                                
+                                if (requirements.requiresUpperBounds()) {
+                                    solverRequirementsData.oneStepTargetProbabilities = computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates);
+                                }
                                 
                                 // Since we needed the transitions to target states to be translated as well for the computation
                                 // of the scheduler, we have to get rid of them now.
-                                eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates);
+                                eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates);
                             }
                         }
                         
@@ -543,7 +564,14 @@ namespace storm {
                         std::vector<ValueType> x(explicitRepresentation.first.getRowGroupCount(), storm::utility::zero<ValueType>());
                         
                         // Now solve the resulting equation system.
-                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create();
+                        
+                        // If the solver requires upper bounds, compute them now.
+                        if (requirements.requiresUpperBounds()) {
+                            setUpperRewardBounds(*solver, dir, explicitRepresentation.first, explicitRepresentation.second, solverRequirementsData.oneStepTargetProbabilities.get());
+                        }
+                        
+                        solver->setMatrix(std::move(explicitRepresentation.first));
                         
                         // Move the scheduler to the solver.
                         if (solverRequirementsData.initialScheduler) {

src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp

@@ -30,6 +30,19 @@ namespace storm {
                 
                 // (3) Compute number of states not in MECs.
                 numberOfMaybeStatesNotInEc = maybeStateToEc.size() - numberOfMaybeStatesInEc;
+                
+                // (4) Compute the number of maybe states that are not in ECs before every maybe state.
+                maybeStatesNotInEcBefore = std::vector<uint64_t>(maybeStateToEc.size());
+                
+                uint64_t count = 0;
+                auto resultIt = maybeStatesNotInEcBefore.begin();
+                for (auto const& e : maybeStateToEc) {
+                    *resultIt = count;
+                    if (e == NOT_IN_EC) {
+                        ++count;
+                    }
+                    ++resultIt;
+                }
             }
 
             template<typename ValueType>
@@ -43,20 +56,8 @@ namespace storm {
             }
             
             template<typename ValueType>
-            std::vector<uint64_t> SparseMdpEndComponentInformation<ValueType>::getNumberOfMaybeStatesNotInEcBeforeIndices() const {
-                std::vector<uint64_t> result(maybeStateToEc.size());
-                
-                uint64_t count = 0;
-                auto resultIt = result.begin();
-                for (auto const& e : maybeStateToEc) {
-                    *resultIt = count;
-                    if (e != NOT_IN_EC) {
-                        ++count;
-                    }
-                    ++resultIt;
-                }
-                
-                return result;
+            std::vector<uint64_t> const& SparseMdpEndComponentInformation<ValueType>::getNumberOfMaybeStatesNotInEcBeforeIndices() const {
+                return maybeStatesNotInEcBefore;
             }
             
             template<typename ValueType>
@@ -71,7 +72,11 @@ namespace storm {
             
             template<typename ValueType>
             uint64_t SparseMdpEndComponentInformation<ValueType>::getRowGroupAfterElimination(uint64_t state) const {
+                if (this->isStateInEc(state)) {
                     return numberOfMaybeStatesNotInEc + getEc(state);
+                } else {
+                    return maybeStatesNotInEcBefore[maybeStatesBefore[state]];
+                }
             }
             
             template<typename ValueType>
@@ -90,20 +95,22 @@ namespace storm {
                 SparseMdpEndComponentInformation<ValueType> result(endComponentDecomposition, maybeStates);
                 
                 // (1) Compute the number of maybe states not in ECs before any other maybe state.
-                std::vector<uint64_t> maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices();
+                std::vector<uint64_t> const& maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices();
                 uint64_t numberOfStates = result.numberOfMaybeStatesNotInEc + result.numberOfEc;
                 
-                STORM_LOG_TRACE("Found " << numberOfStates << " states, " << result.numberOfMaybeStatesNotInEc << " not in ECs, " << result.numberOfMaybeStatesInEc << " in ECs and " << result.numberOfEc << " ECs.");
+                STORM_LOG_TRACE("Creating " << numberOfStates << " states from " << result.numberOfMaybeStatesNotInEc << " states not in ECs and " << result.numberOfMaybeStatesInEc << " states in " << result.numberOfEc << " ECs.");
                 
                 // Prepare builder and vector storage.
                 storm::storage::SparseMatrixBuilder<ValueType> builder(0, numberOfStates, 0, true, true, numberOfStates);
                 STORM_LOG_ASSERT((sumColumns && columnSumVector) || (!sumColumns && !columnSumVector), "Expecting a bit vector for which columns to sum iff there is a column sum result vector.");
                 if (columnSumVector) {
-                    columnSumVector->resize(numberOfStates);
+                    // Clearing column sum vector as we do not know the number of choices at this point.
+                    columnSumVector->resize(0);
                 }
                 STORM_LOG_ASSERT((summand && summandResultVector) || (!summand && !summandResultVector), "Expecting summand iff there is a summand result vector.");
                 if (summandResultVector) {
-                    summandResultVector->resize(numberOfStates);
+                    // Clearing summand result vector as we do not know the number of choices at this point.
+                    summandResultVector->resize(0);
                 }
                 std::vector<std::pair<uint64_t, ValueType>> ecValuePairs;
                 
@@ -121,14 +128,17 @@ namespace storm {
                             ecValuePairs.clear();
                             
                             if (summand) {
-                                (*summandResultVector)[currentRow] += (*summand)[row];
+                                summandResultVector->emplace_back((*summand)[row]);
+                            }
+                            if (columnSumVector) {
+                                columnSumVector->emplace_back(storm::utility::zero<ValueType>());
                             }
                             for (auto const& e : transitionMatrix.getRow(row)) {
                                 if (sumColumns && sumColumns->get(e.getColumn())) {
-                                    (*columnSumVector)[currentRow] += e.getValue();
+                                    columnSumVector->back() += e.getValue();
                                 } else if (maybeStates.get(e.getColumn())) {
                                     // If the target state of the transition is not contained in an EC, we can just add the entry.
-                                    if (result.isStateInEc(e.getColumn())) {
+                                    if (!result.isStateInEc(e.getColumn())) {
                                         builder.addNextValue(currentRow, maybeStatesNotInEcBefore[result.maybeStatesBefore[e.getColumn()]], e.getValue());
                                     } else {
                                         // Otherwise, we store the information that the state can go to a certain EC.
@@ -170,14 +180,17 @@ namespace storm {
                             ecValuePairs.clear();
                             
                             if (summand) {
-                                (*summandResultVector)[currentRow] += (*summand)[row];
+                                summandResultVector->emplace_back((*summand)[row]);
+                            }
+                            if (columnSumVector) {
+                                columnSumVector->emplace_back(storm::utility::zero<ValueType>());
                             }
                             for (auto const& e : transitionMatrix.getRow(row)) {
                                 if (sumColumns && sumColumns->get(e.getColumn())) {
-                                    (*columnSumVector)[currentRow] += e.getValue();
+                                    columnSumVector->back() += e.getValue();
                                 } else if (maybeStates.get(e.getColumn())) {
                                     // If the target state of the transition is not contained in an EC, we can just add the entry.
-                                    if (result.isStateInEc(e.getColumn())) {
+                                    if (!result.isStateInEc(e.getColumn())) {
                                         builder.addNextValue(currentRow, maybeStatesNotInEcBefore[result.maybeStatesBefore[e.getColumn()]], e.getValue());
                                     } else {
                                         // Otherwise, we store the information that the state can go to a certain EC.

src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h

@@ -29,7 +29,7 @@ namespace storm {
                  * Retrieves for each maybe state the number of maybe states not contained in ECs with an index smaller
                  * than the requested one.
                  */
-                std::vector<uint64_t> getNumberOfMaybeStatesNotInEcBeforeIndices() const;
+                std::vector<uint64_t> const& getNumberOfMaybeStatesNotInEcBeforeIndices() const;
                 
                 /*!
                  * Retrieves the total number of maybe states on in ECs.
@@ -65,6 +65,7 @@ namespace storm {
 
                 // Data about end components.
                 std::vector<uint64_t> maybeStatesBefore;
+                std::vector<uint64_t> maybeStatesNotInEcBefore;
                 uint64_t numberOfMaybeStatesInEc;
                 uint64_t numberOfMaybeStatesNotInEc;
                 uint64_t numberOfEc;

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -449,7 +449,7 @@ namespace storm {
             }
             
             template<typename ValueType>
-            void computeFixedPointSystemUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
+            void computeFixedPointSystemUntilProbabilities(storm::solver::SolveGoal<ValueType>& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
                 // First, we can eliminate the rows and columns from the original transition probability matrix for states
                 // whose probabilities are already known.
                 submatrix = transitionMatrix.getSubmatrix(true, qualitativeStateSets.maybeStates, qualitativeStateSets.maybeStates, false);
@@ -457,10 +457,13 @@ namespace storm {
                 // Prepare the right-hand side of the equation system. For entry i this corresponds to
                 // the accumulated probability of going from state i to some state that has probability 1.
                 b = transitionMatrix.getConstrainedRowGroupSumVector(qualitativeStateSets.maybeStates, qualitativeStateSets.statesWithProbability1);
+                
+                // If the solve goal has relevant values, we need to adjust them.
+                goal.restrictRelevantValues(qualitativeStateSets.maybeStates);
             }
             
             template<typename ValueType>
-            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
+            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(storm::solver::SolveGoal<ValueType>& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b) {
                 
                 // Start by computing the states that are in MECs.
                 storm::storage::MaximalEndComponentDecomposition<ValueType> endComponentDecomposition(transitionMatrix, backwardTransitions, qualitativeStateSets.maybeStates);
@@ -468,10 +471,26 @@ namespace storm {
                 // Only do more work if there are actually end-components.
                 if (!endComponentDecomposition.empty()) {
                     STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " EC(s).");
-                    return SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, &b, nullptr);
+                    SparseMdpEndComponentInformation<ValueType> result = SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, &b, nullptr);
+                    
+                    // If the solve goal has relevant values, we need to adjust them.
+                    if (goal.hasRelevantValues()) {
+                        storm::storage::BitVector newRelevantValues(submatrix.getRowGroupCount());
+                        for (auto state : goal.relevantValues()) {
+                            if (qualitativeStateSets.maybeStates.get(state)) {
+                                newRelevantValues.set(result.getRowGroupAfterElimination(state));
+                            }
+                        }
+                        if (!newRelevantValues.empty()) {
+                            goal.setRelevantValues(std::move(newRelevantValues));
+                        }
+                    }
+                    
+                    return result;
                 } else {
                     STORM_LOG_DEBUG("Not eliminating ECs as there are none.");
-                    computeFixedPointSystemUntilProbabilities(transitionMatrix, qualitativeStateSets, submatrix, b);
+                    computeFixedPointSystemUntilProbabilities(goal, transitionMatrix, qualitativeStateSets, submatrix, b);
+                    
                     return boost::none;
                 }
             }
@@ -516,17 +535,16 @@ namespace storm {
                         // If the hint information tells us that we have to eliminate MECs, we do so now.
                         boost::optional<SparseMdpEndComponentInformation<ValueType>> ecInformation;
                         if (hintInformation.getEliminateEndComponents()) {
-                            ecInformation = computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, submatrix, b);
+                            ecInformation = computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(goal, transitionMatrix, backwardTransitions, qualitativeStateSets, submatrix, b);
 
                             // Make sure we are not supposed to produce a scheduler if we actually eliminate end components.
                             STORM_LOG_THROW(!ecInformation || !ecInformation.get().getEliminatedEndComponents() || !produceScheduler, storm::exceptions::NotSupportedException, "Producing schedulers is not supported if end-components need to be eliminated for the solver.");
                         } else {
                             // Otherwise, we compute the standard equations.
-                            computeFixedPointSystemUntilProbabilities(transitionMatrix, qualitativeStateSets, submatrix, b);
+                            computeFixedPointSystemUntilProbabilities(goal, transitionMatrix, qualitativeStateSets, submatrix, b);
                         }
                         
                         // Now compute the results for the maybe states.
-                        goal.restrictRelevantValues(qualitativeStateSets.maybeStates);
                         MaybeStateResult<ValueType> resultForMaybeStates = computeValuesForMaybeStates(std::move(goal), std::move(submatrix), b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation);
                         
                         // If we eliminated end components, we need to extract the result differently.
@@ -725,7 +743,7 @@ namespace storm {
             }
             
             template<typename ValueType>
-            void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, std::vector<ValueType>* oneStepTargetProbabilities = nullptr) {
+            void computeFixedPointSystemReachabilityRewards(storm::solver::SolveGoal<ValueType>& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, std::vector<ValueType>* oneStepTargetProbabilities = nullptr) {
                 // Remove rows and columns from the original transition probability matrix for states whose reward values are already known.
                 // If there are infinity states, we additionally have to remove choices of maybeState that lead to infinity.
                 if (qualitativeStateSets.infinityStates.empty()) {
@@ -742,10 +760,13 @@ namespace storm {
                         (*oneStepTargetProbabilities) = transitionMatrix.getConstrainedRowSumVector(*selectedChoices, targetStates);
                     }
                 }
+                
+                // If the solve goal has relevant values, we need to adjust them.
+                goal.restrictRelevantValues(qualitativeStateSets.maybeStates);
             }
             
             template<typename ValueType>
-            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, boost::optional<std::vector<ValueType>>& oneStepTargetProbabilities) {
+            boost::optional<SparseMdpEndComponentInformation<ValueType>> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::solver::SolveGoal<ValueType>& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional<storm::storage::BitVector> const& selectedChoices, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix<ValueType>& submatrix, std::vector<ValueType>& b, boost::optional<std::vector<ValueType>>& oneStepTargetProbabilities) {
                 
                 // Start by computing the choices with reward 0, as we only want ECs within this fragment.
                 storm::storage::BitVector zeroRewardChoices(transitionMatrix.getRowCount());
@@ -789,10 +810,25 @@ namespace storm {
                 // Only do more work if there are actually end-components.
                 if (doDecomposition && !endComponentDecomposition.empty()) {
                     STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " ECs.");
-                    return SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, oneStepTargetProbabilities ? &targetStates : nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr, &b);
+                    SparseMdpEndComponentInformation<ValueType> result = SparseMdpEndComponentInformation<ValueType>::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, oneStepTargetProbabilities ? &targetStates : nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr, &b);
+                    
+                    // If the solve goal has relevant values, we need to adjust them.
+                    if (goal.hasRelevantValues()) {
+                        storm::storage::BitVector newRelevantValues(submatrix.getRowGroupCount());
+                        for (auto state : goal.relevantValues()) {
+                            if (qualitativeStateSets.maybeStates.get(state)) {
+                                newRelevantValues.set(result.getRowGroupAfterElimination(state));
+                            }
+                        }
+                        if (!newRelevantValues.empty()) {
+                            goal.setRelevantValues(std::move(newRelevantValues));
+                        }
+                    }
+                    
+                    return result;
                 } else {
                     STORM_LOG_DEBUG("Not eliminating ECs as there are none.");
-                    computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr);
+                    computeFixedPointSystemReachabilityRewards(goal, transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr);
                     return boost::none;
                 }
             }
@@ -862,13 +898,13 @@ namespace storm {
                         // If the hint information tells us that we have to eliminate MECs, we do so now.
                         boost::optional<SparseMdpEndComponentInformation<ValueType>> ecInformation;
                         if (hintInformation.getEliminateEndComponents()) {
-                            ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities);
+                            ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(goal, transitionMatrix, backwardTransitions, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities);
                             
                             // Make sure we are not supposed to produce a scheduler if we actually eliminate end components.
                             STORM_LOG_THROW(!ecInformation || !ecInformation.get().getEliminatedEndComponents() || !produceScheduler, storm::exceptions::NotSupportedException, "Producing schedulers is not supported if end-components need to be eliminated for the solver.");
                         } else {
                             // Otherwise, we compute the standard equations.
-                            computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr);
+                            computeFixedPointSystemReachabilityRewards(goal, transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr);
                         }
                         
                         // If we need to compute upper bounds, do so now.
@@ -878,7 +914,6 @@ namespace storm {
                         }
                         
                         // Now compute the results for the maybe states.
-                        goal.restrictRelevantValues(qualitativeStateSets.maybeStates);
                         MaybeStateResult<ValueType> resultForMaybeStates = computeValuesForMaybeStates(std::move(goal), std::move(submatrix), b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation);
 
                         // If we eliminated end components, we need to extract the result differently.

src/storm/models/symbolic/StandardRewardModel.cpp

@@ -101,6 +101,21 @@ namespace storm {
                 return result;
             }
             
+            template <storm::dd::DdType Type, typename ValueType>
+            storm::dd::Add<Type, ValueType> StandardRewardModel<Type, ValueType>::getTotalRewardVector(storm::dd::Add<Type, ValueType> const& stateFilterAdd, storm::dd::Add<Type, ValueType> const& choiceFilterAdd, storm::dd::Add<Type, ValueType> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const {
+                storm::dd::Add<Type, ValueType> result = transitionMatrix.getDdManager().template getAddZero<ValueType>();
+                if (this->hasStateRewards()) {
+                    result += stateFilterAdd * choiceFilterAdd * optionalStateRewardVector.get();
+                }
+                if (this->hasStateActionRewards()) {
+                    result += choiceFilterAdd * optionalStateActionRewardVector.get();
+                }
+                if (this->hasTransitionRewards()) {
+                    result += (transitionMatrix * this->getTransitionRewardMatrix()).sumAbstract(columnVariables);
+                }
+                return result;
+            }
+            
             template <storm::dd::DdType Type, typename ValueType>
             storm::dd::Add<Type, ValueType> StandardRewardModel<Type, ValueType>::getTotalRewardVector(storm::dd::Add<Type, ValueType> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const {
                 storm::dd::Add<Type, ValueType> result = transitionMatrix.getDdManager().template getAddZero<ValueType>();

src/storm/models/symbolic/StandardRewardModel.h

@@ -161,6 +161,19 @@ namespace storm {
                  */
                 storm::dd::Add<Type, ValueType> getTotalRewardVector(storm::dd::Add<Type, ValueType> const& filterAdd, storm::dd::Add<Type, ValueType> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const;
                 
+                /*!
+                 * Creates a vector representing the complete reward vector based on the state-, state-action- and
+                 * transition-based rewards in the reward model. The state- and state-action rewards are filtered with
+                 * the given filter DD.
+                 *
+                 * @param stateFilterAdd The DD used for filtering the state rewards.
+                 * @param choiceFilterAdd The DD used for filtering the state action rewards.
+                 * @param transitionMatrix The matrix that is used to weight the values of the transition reward matrix.
+                 * @param columnVariables The column variables of the model.
+                 * @return The full state-action reward vector.
+                 */
+                storm::dd::Add<Type, ValueType> getTotalRewardVector(storm::dd::Add<Type, ValueType> const& stateFilterAdd, storm::dd::Add<Type, ValueType> const& choiceFilterAdd, storm::dd::Add<Type, ValueType> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const;
+                
                 /*!
                  * Creates a vector representing the complete reward vector based on the state-, state-action- and
                  * transition-based rewards in the reward model.

src/storm/solver/GmmxxMultiplier.cpp

@@ -80,11 +80,19 @@ namespace storm {
             
             uint64_t choice;
             for (auto row_group_it = rowGroupIndices.end() - 2, row_group_ite = rowGroupIndices.begin() - 1; row_group_it != row_group_ite; --row_group_it, --choice_it, --target_it) {
-                T currentValue = b ? *add_it : storm::utility::zero<T>();
-                --add_it;
+                if (choices) {
+                    *choice_it = 0;
+                }
+
+                T currentValue = storm::utility::zero<T>();
                 
                 // Only multiply and reduce if the row group is not empty.
                 if (*row_group_it != *(row_group_it + 1)) {
+                    if (b) {
+                        currentValue = *add_it;
+                        --add_it;
+                    }
+                    
                     currentValue += vect_sp(gmm::linalg_traits<MatrixType>::row(itr), x);
                     
                     if (choices) {
@@ -94,7 +102,7 @@ namespace storm {
                     
                     --itr;
                     
-                    for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr) {
+                    for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr, --add_it) {
                         T newValue = b ? *add_it : storm::utility::zero<T>();
                         newValue += vect_sp(gmm::linalg_traits<MatrixType>::row(itr), x);
                         
@@ -108,9 +116,6 @@ namespace storm {
                                 *choice_it = choice;
                             }
                         }
-                        if (b) {
-                            --add_it;
-                        }
                     }
                 } else if (choices) {
                     *choice_it = 0;
@@ -165,14 +170,19 @@ namespace storm {
                 auto resultIt = result.begin() + range.begin();
 
                 for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) {
-                    T currentValue = b ? *bIt : storm::utility::zero<T>();
-                    ++bIt;
                     if (choices) {
                         *choiceIt = 0;
                     }
                     
+                    T currentValue = storm::utility::zero<T>();
+                    
                     // Only multiply and reduce if the row group is not empty.
                     if (*groupIt != *(groupIt + 1)) {
+                        if (b) {
+                            currentValue = *bIt;
+                            ++bIt;
+                        }
+                        
                         ++itr;
                         
                         for (auto itre = mat_row_const_begin(matrix) + *(groupIt + 1); itr != itre; ++itr) {

src/storm/solver/SolveGoal.cpp

@@ -117,6 +117,11 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void SolveGoal<ValueType>::setRelevantValues(storm::storage::BitVector&& values) {
+            relevantValueVector = std::move(values);
+        }
+
         template class SolveGoal<double>;
         
 #ifdef STORM_HAVE_CARL

src/storm/solver/SolveGoal.h

@@ -82,6 +82,7 @@ namespace storm {
             storm::storage::BitVector const& relevantValues() const;
            
             void restrictRelevantValues(storm::storage::BitVector const& filter);
+            void setRelevantValues(storm::storage::BitVector&& values);
             
         private:
             boost::optional<OptimizationDirection> optimizationDirection;

src/storm/storage/SparseMatrix.cpp

@@ -1565,20 +1565,21 @@ namespace storm {
             }
             
             for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) {
-                ValueType currentValue = summand ? *summandIt : storm::utility::zero<ValueType>();
-                ++summandIt;
+                ValueType currentValue = storm::utility::zero<ValueType>();
                 if (choices) {
                     *choiceIt = 0;
                 }
                 
                 // 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) {
-                        *choiceIt = 0;
-                    }
                     
                     ++rowIt;
                     
@@ -1598,6 +1599,8 @@ namespace storm {
                             ++summandIt;
                         }
                     }
+                } else if (choices) {
+                    *choiceIt = 0;
                 }
                 
                 // Finally write value to target vector.
@@ -1627,21 +1630,27 @@ namespace storm {
             }
             
             for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt) {
-                ValueType currentValue = summand ? *summandIt : storm::utility::zero<ValueType>();
-                --summandIt;
+                if (choices) {
+                    *choiceIt = 0;
+                }
+                ValueType currentValue = storm::utility::zero<ValueType>();
                 
                 // 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) {
                         *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt;
                     }
-                    
                     --rowIt;
                     
-                    for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, ++i) {
+                    for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, ++i, --summandIt) {
                         ValueType newValue = summand ? *summandIt : storm::utility::zero<ValueType>();
                         for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) {
                             newValue += elementIt->getValue() * vector[elementIt->getColumn()];
@@ -1653,13 +1662,8 @@ namespace storm {
                                 *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt;
                             }
                         }
-                        if (summand) {
-                            --summandIt;
                     }
                 }
-                } else if (choices) {
-                    *choiceIt = 0;
-                }
                 
                 // Finally write value to target vector.
                 *resultIt = currentValue;
@@ -1703,14 +1707,19 @@ namespace storm {
                 auto resultIt = result.begin() + range.begin();
                 
                 for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) {
-                    ValueType currentValue = summand ? *summandIt : storm::utility::zero<ValueType>();
-                    ++summandIt;
                     if (choices) {
                         *choiceIt = 0;
                     }
                     
+                    ValueType currentValue = storm::utility::zero<ValueType>();
+                    
                     // Only multiply and reduce if there is at least one row in the group.
                     if (*groupIt < *(groupIt + 1)) {
+                        if (summand) {
+                            currentValue = *summandIt;
+                            ++summandIt;
+                        }
+
                         for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) {
                             currentValue += elementIt->getValue() * x[elementIt->getColumn()];
                         }

src/storm/storage/dd/Add.cpp

@@ -473,6 +473,67 @@ namespace storm {
             return result;
         }
         
+        template<DdType LibraryType, typename ValueType>
+        std::vector<ValueType> Add<LibraryType, ValueType>::toVector(storm::dd::Add<LibraryType, ValueType> const& matrix, std::vector<uint_fast64_t> const& rowGroupIndices, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const {
+            std::vector<uint_fast64_t> ddRowVariableIndices;
+            std::vector<uint_fast64_t> ddColumnVariableIndices;
+            std::vector<uint_fast64_t> ddGroupVariableIndices;
+            std::set<storm::expressions::Variable> rowAndColumnMetaVariables;
+            
+            for (auto const& variable : rowMetaVariables) {
+                DdMetaVariable<LibraryType> const& metaVariable = this->getDdManager().getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddRowVariableIndices.push_back(ddVariable.getIndex());
+                }
+                rowAndColumnMetaVariables.insert(variable);
+            }
+            std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
+            for (auto const& variable : columnMetaVariables) {
+                DdMetaVariable<LibraryType> const& metaVariable = this->getDdManager().getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddColumnVariableIndices.push_back(ddVariable.getIndex());
+                }
+                rowAndColumnMetaVariables.insert(variable);
+            }
+            std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
+            for (auto const& variable : groupMetaVariables) {
+                DdMetaVariable<LibraryType> const& metaVariable = this->getDdManager().getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddGroupVariableIndices.push_back(ddVariable.getIndex());
+                }
+            }
+            std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end());
+            
+            Bdd<LibraryType> columnVariableCube = Bdd<LibraryType>::getCube(this->getDdManager(), columnMetaVariables);
+            
+            // Copy the row group indices so we can modify them.
+            std::vector<uint_fast64_t> mutableRowGroupIndices = rowGroupIndices;
+            
+            // Create the explicit vector we need to fill later.
+            std::vector<ValueType> explicitVector(mutableRowGroupIndices.back());
+            
+            // Next, we split the matrix into one for each group. Note that this only works if the group variables are at the very top.
+            std::vector<std::pair<InternalAdd<LibraryType, ValueType>, InternalAdd<LibraryType, ValueType>>> internalAddGroups = matrix.internalAdd.splitIntoGroups(*this, ddGroupVariableIndices);
+            std::vector<std::pair<Add<LibraryType, ValueType>, Add<LibraryType, ValueType>>> groups;
+            for (auto const& internalAdd : internalAddGroups) {
+                groups.push_back(std::make_pair(Add<LibraryType, ValueType>(this->getDdManager(), internalAdd.first, rowAndColumnMetaVariables), Add<LibraryType, ValueType>(this->getDdManager(), internalAdd.second, rowMetaVariables)));
+            }
+            
+            std::vector<InternalAdd<LibraryType, uint_fast64_t>> statesWithGroupEnabled(groups.size());
+            for (uint_fast64_t i = 0; i < groups.size(); ++i) {
+                std::pair<Add<LibraryType, ValueType>, Add<LibraryType, ValueType>> const& ddPair = groups[i];
+                Bdd<LibraryType> matrixDdNotZero = ddPair.first.notZero();
+                Bdd<LibraryType> vectorDdNotZero = ddPair.second.notZero();
+                
+                ddPair.second.internalAdd.composeWithExplicitVector(rowOdd, ddRowVariableIndices, mutableRowGroupIndices, explicitVector, std::plus<ValueType>());
+                
+                InternalAdd<LibraryType, uint_fast64_t> statesWithGroupEnabled = (matrixDdNotZero.existsAbstract(columnMetaVariables) || vectorDdNotZero).template toAdd<uint_fast64_t>();
+                statesWithGroupEnabled.composeWithExplicitVector(rowOdd, ddRowVariableIndices, mutableRowGroupIndices, std::plus<uint_fast64_t>());
+            }
+            
+            return explicitVector;
+        }
+                
         template<DdType LibraryType, typename ValueType>
         storm::storage::SparseMatrix<ValueType> Add<LibraryType, ValueType>::toMatrix() const {
             std::set<storm::expressions::Variable> rowVariables;
@@ -697,7 +758,7 @@ namespace storm {
         }
 
         template<DdType LibraryType, typename ValueType>
-        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> Add<LibraryType, ValueType>::toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const {
+        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> Add<LibraryType, ValueType>::toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const {
             std::set<storm::expressions::Variable> rowMetaVariables;
             std::set<storm::expressions::Variable> columnMetaVariables;
             
@@ -715,11 +776,11 @@ namespace storm {
             }
             
             // Create the canonical row group sizes and build the matrix.
-            return toMatrixVector(vector, std::move(rowGroupSizes), rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd);
+            return toMatrixVector(vector, rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd);
         }
 
         template<DdType LibraryType, typename ValueType>
-        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> Add<LibraryType, ValueType>::toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::vector<uint_fast64_t>&& rowGroupIndices, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const {
+        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> Add<LibraryType, ValueType>::toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const {
             std::vector<uint_fast64_t> ddRowVariableIndices;
             std::vector<uint_fast64_t> ddColumnVariableIndices;
             std::vector<uint_fast64_t> ddGroupVariableIndices;
@@ -751,6 +812,9 @@ namespace storm {
             
             Bdd<LibraryType> columnVariableCube = Bdd<LibraryType>::getCube(this->getDdManager(), columnMetaVariables);
 
+            // Count how many choices each row group has.
+            std::vector<uint_fast64_t> rowGroupIndices = (this->notZero().existsAbstract(columnMetaVariables) || vector.notZero()).template toAdd<uint_fast64_t>().sumAbstract(groupMetaVariables).toVector(rowOdd);
+            
             // Transform the row group sizes to the actual row group indices.
             rowGroupIndices.resize(rowGroupIndices.size() + 1);
             uint_fast64_t tmp = 0;

src/storm/storage/dd/Add.h

@@ -556,7 +556,26 @@ namespace storm {
             std::vector<ValueType> toVector(storm::dd::Odd const& rowOdd) const;
             
             /*!
-             * Converts the ADD to a (sparse) double matrix. All contained non-primed variables are assumed to encode the
+             * Converts the ADD to a row-grouped vector while respecting the row group sizes of the provided matrix.
+             * That is, if the vector has a zero entry for some row in a row group for which the matrix has a non-zero
+             * row, the value at the vector will be correctly set to zero. Note: this function assumes that the meta
+             * variables used to distinguish different row groups are at the very top of the ADD.
+             *
+             * @param matrix The symbolic matrix whose row group sizes to respect.
+             * @param rowGroupSizes A vector specifying the sizes of the row groups.
+             * @param rowMetaVariables The meta variables that encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables that encode the columns of the matrix.
+             * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
+             * @param rowOdd The ODD used for determining the correct row.
+             * @param columnOdd The ODD used for determining the correct column.
+             * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector
+             * (if it was given).
+             * @return The vector that is represented by this ADD.
+             */
+            std::vector<ValueType> toVector(storm::dd::Add<LibraryType, ValueType> const& matrix, std::vector<uint_fast64_t> const& rowGroupSizes, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
+            
+            /*!
+             * Converts the ADD to a (sparse) matrix. All contained non-primed variables are assumed to encode the
              * row, whereas all primed variables are assumed to encode the column.
              *
              * @return The matrix that is represented by this ADD.
@@ -564,7 +583,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> toMatrix() const;
             
             /*!
-             * Converts the ADD to a (sparse) double matrix. All contained non-primed variables are assumed to encode the
+             * Converts the ADD to a (sparse) matrix. All contained non-primed variables are assumed to encode the
              * row, whereas all primed variables are assumed to encode the column. The given offset-labeled DDs are used
              * to determine the correct row and column, respectively, for each entry.
              *
@@ -575,7 +594,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> toMatrix(storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
             
             /*!
-             * Converts the ADD to a (sparse) double matrix. The given offset-labeled DDs are used to determine the
+             * Converts the ADD to a (sparse) matrix. The given offset-labeled DDs are used to determine the
              * correct row and column, respectively, for each entry.
              *
              * @param rowMetaVariables The meta variables that encode the rows of the matrix.
@@ -587,7 +606,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> toMatrix(std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
             
             /*!
-             * Converts the ADD to a row-grouped (sparse) double matrix. The given offset-labeled DDs are used to
+             * Converts the ADD to a row-grouped (sparse) matrix. The given offset-labeled DDs are used to
              * determine the correct row and column, respectively, for each entry. Note: this function assumes that
              * the meta variables used to distinguish different row groups are at the very top of the ADD.
              *
@@ -599,19 +618,18 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> toMatrix(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
             
             /*!
-             * Converts the ADD to a row-grouped (sparse) double matrix and the given vector to a row-grouped vector.
+             * Converts the ADD to a row-grouped (sparse) matrix and the given vector to a row-grouped vector.
              * The given offset-labeled DDs are used to determine the correct row and column, respectively, for each
              * entry. Note: this function assumes that the meta variables used to distinguish different row groups are
              * at the very top of the ADD.
              *
              * @param vector The symbolic vector to convert.
-             * @param rowGroupSizes A vector specifying the sizes of the row groups.
              * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
              * @param rowOdd The ODD used for determining the correct row.
              * @param columnOdd The ODD used for determining the correct column.
              * @return The matrix that is represented by this ADD.
              */
-            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
+            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
             
             /*!
              * Exports the DD to the given file in the dot format.
@@ -702,7 +720,6 @@ namespace storm {
              * different row groups are at the very top of the ADD.
              *
              * @param vector The vector that is to be transformed to an equally grouped explicit vector.
-             * @param rowGroupSizes A vector specifying the sizes of the row groups.
              * @param rowMetaVariables The meta variables that encode the rows of the matrix.
              * @param columnMetaVariables The meta variables that encode the columns of the matrix.
              * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
@@ -711,7 +728,7 @@ namespace storm {
              * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector
              * (if it was given).
              */
-            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
+            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> toMatrixVector(storm::dd::Add<LibraryType, ValueType> const& vector, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const;
                         
             // The internal ADD that depends on the chosen library.
             InternalAdd<LibraryType, ValueType> internalAdd;

src/storm/utility/vector.h

@@ -622,7 +622,9 @@ namespace storm {
                         choiceIt = choices->begin() + startRow;
                     }
                     
-                    for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) {
+                    for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) {
+                        // Only traverse elements if the row group is non-empty.
+                        if (*rowGroupingIt != *(rowGroupingIt + 1)) {
                             *targetIt = *sourceIt;
                             ++sourceIt;
                             localChoice = 1;
@@ -638,9 +640,8 @@ namespace storm {
                                     }
                                 }
                             }
-                        
-                        if (choices != nullptr) {
-                            ++choiceIt;
+                        } else {
+                            *targetIt = storm::utility::zero<T>();
                         }
                     }
                 }
@@ -678,7 +679,9 @@ namespace storm {
                     choiceIt = choices->begin();
                 }
                 
-                for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) {
+                for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) {
+                    // Only traverse elements if the row group is non-empty.
+                    if (*rowGroupingIt != *(rowGroupingIt + 1)) {
                         *targetIt = *sourceIt;
                         ++sourceIt;
                         localChoice = 1;
@@ -693,8 +696,8 @@ namespace storm {
                                 }
                             }
                         }
-                    if (choices != nullptr) {
-                        ++choiceIt;
+                    } else {
+                        *targetIt = storm::utility::zero<T>();
                     }
                 }
             }