moved multidimensional unfolding implementation from multiobjective into helper namespace

7 years ago · 86253fe88a
15 changed files with 1999 additions and 17 deletions
--- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h
@ -4,7 +4,7 @@
 #include <vector>

 #include "storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h"
-#include "storm/modelchecker/multiobjective/rewardbounded/MultiDimensionalRewardUnfolding.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"

 namespace storm {
    namespace modelchecker {
--- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp
@ -84,7 +84,6 @@ namespace storm {
                    storm::utility::exportDataToCSVFile<ValueType, ValueType, std::string>("cdf" + std::to_string(numChecks) + ".csv", cdfData, weightVector, headers);
                }
                auto solution = rewardUnfolding.getInitialStateResult(initEpoch);
-                // Todo: we currently assume precise results...
                auto solutionIt = solution.begin();
                ++solutionIt;
                underApproxResult = std::vector<ValueType>(solutionIt, solution.end());
@ -93,14 +92,14 @@ namespace storm {
            }
            
            template <class SparseMdpModelType>
-            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision) {
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::computeEpochSolution(typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision) {
                
                ++numCheckedEpochs;
                swEpochModelBuild.start();
                auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
                swEpochModelBuild.stop();
                swEpochModelAnalysis.start();
-                std::vector<typename MultiDimensionalRewardUnfolding<ValueType, false>::SolutionType> result;
+                std::vector<typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::SolutionType> result;
                result.reserve(epochModel.epochInStates.getNumberOfSetBits());
                uint64_t solutionSize = this->objectives.size() + 1;

@ -287,7 +286,7 @@ namespace storm {
            }

            template <class SparseMdpModelType>
-            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision) {
+            void SparseMdpRewardBoundedPcaaWeightVectorChecker<SparseMdpModelType>::updateCachedData(typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision) {
                if (epochModel.epochMatrixChanged) {
                
                    // Update the cached MinMaxSolver data
--- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h
@ -3,7 +3,7 @@
 #include <vector>

 #include "storm/modelchecker/multiobjective/pcaa/PcaaWeightVectorChecker.h"
-#include "storm/modelchecker/multiobjective/rewardbounded/MultiDimensionalRewardUnfolding.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
 #include "storm/solver/MinMaxLinearEquationSolver.h"
 #include "storm/solver/LinearEquationSolver.h"

@ -70,17 +70,17 @@ namespace storm {
                    std::vector<std::vector<ValueType>> xLinEq;
                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSolver;
                    
-                    std::vector<typename MultiDimensionalRewardUnfolding<ValueType, false>::SolutionType> solutions;
+                    std::vector<typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::SolutionType> solutions;
                };
                
-                void computeEpochSolution(typename MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision);
+                void computeEpochSolution(typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::Epoch const& epoch, std::vector<ValueType> const& weightVector, EpochCheckingData& cachedData, ValueType const& precision);
                
-                void updateCachedData(typename MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision);
+                void updateCachedData(typename helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false>::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector<ValueType> const& weightVector, ValueType const& precision);
                
                storm::utility::Stopwatch swAll, swEpochModelBuild, swEpochModelAnalysis;
                uint64_t numSchedChanges, numCheckedEpochs, numChecks;
                
-                MultiDimensionalRewardUnfolding<ValueType, false> rewardUnfolding;
+                helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, false> rewardUnfolding;
                
                boost::optional<std::vector<ValueType>> underApproxResult;
                boost::optional<std::vector<ValueType>> overApproxResult;
--- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@ -66,7 +66,7 @@ namespace storm {
                    opInfo.bound = checkTask.getBound();
                }
                auto formula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(checkTask.getFormula().asSharedPointer(), opInfo);
-                storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
+                helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(helper::SolutionType::UntilProbabilities, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
                return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
            } else {
@ -145,7 +145,7 @@ namespace storm {
                    opInfo.bound = checkTask.getBound();
                }
                auto formula = std::make_shared<storm::logic::RewardOperatorFormula>(checkTask.getFormula().asSharedPointer(), checkTask.getRewardModel(), opInfo);
-                storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
+                helper::rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(this->getModel(), formula);
                auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(helper::SolutionType::ExpectedRewards, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory);
                return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
            } else {
--- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
@ -84,7 +84,7 @@ namespace storm {

            
            template<typename ValueType>
-            std::vector<ValueType> analyzeTrivialEpochModel(OptimizationDirection dir, typename storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel) {
+            std::vector<ValueType> analyzeTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel) {
                // Assert that the epoch model is indeed trivial
                assert(epochModel.epochMatrix.getEntryCount() == 0);
                
@ -134,7 +134,7 @@ namespace storm {
            }
            
            template<typename ValueType>
-            std::vector<ValueType> analyzeNonTrivialEpochModel(OptimizationDirection dir, typename storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ValueType const& precision, SolutionType const& type) {
+            std::vector<ValueType> analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>::EpochModel& epochModel, std::vector<ValueType>& x, std::vector<ValueType>& b, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ValueType const& precision, SolutionType const& type) {
 
                // Update some data for the case that the Matrix has changed
                if (epochModel.epochMatrixChanged) {
@ -182,7 +182,7 @@ namespace storm {
            }
            
            template<typename ValueType>
-            std::map<storm::storage::sparse::state_type, ValueType> SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+            std::map<storm::storage::sparse::state_type, ValueType> SparseMdpPrctlHelper<ValueType>::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
                storm::utility::Stopwatch swAll(true), swBuild, swCheck;
                auto initEpoch = rewardUnfolding.getStartEpoch();
                auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
--- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h
+++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h
@ -7,7 +7,7 @@
 #include "storm/modelchecker/prctl/helper/SolutionType.h"
 #include "storm/storage/SparseMatrix.h"
 #include "storm/storage/MaximalEndComponent.h"
-#include "storm/modelchecker/multiobjective/rewardbounded/MultiDimensionalRewardUnfolding.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
 #include "MDPModelCheckingHelperReturnType.h"

 #include "storm/utility/solver.h"
@ -38,7 +38,7 @@ namespace storm {
                
                static std::vector<ValueType> computeStepBoundedUntilProbabilities(storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint());

-                static std::map<storm::storage::sparse::state_type, ValueType> computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                static std::map<storm::storage::sparse::state_type, ValueType> computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
                
                static std::vector<ValueType> computeNextProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);

--- a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h
@ -0,0 +1,25 @@
+#pragma once
+
+#include <boost/optional.hpp>
+
+#include "storm/storage/BitVector.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                template<typename ValueType>
+                struct Dimension {
+                    std::shared_ptr<storm::logic::Formula const> formula;
+                    uint64_t objectiveIndex;
+                    boost::optional<std::string> memoryLabel;
+                    bool isUpperBounded;
+                    ValueType scalingFactor;
+                    storm::storage::BitVector dependentDimensions;
+                    boost::optional<uint64_t> maxValue;
+                };
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
@ -0,0 +1,303 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h"
+
+#include "storm/utility/macros.h"
+
+#include "storm/exceptions/IllegalArgumentException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                EpochManager::EpochManager() : dimensionCount(0) {
+                    // Intentionally left empty
+                }
+                
+                EpochManager::EpochManager(uint64_t dimensionCount) : dimensionCount(dimensionCount) {
+                    STORM_LOG_THROW(dimensionCount > 0, storm::exceptions::IllegalArgumentException, "Invoked EpochManager with zero dimension count.");
+                    STORM_LOG_THROW(dimensionCount <= 64, storm::exceptions::IllegalArgumentException, "Invoked EpochManager with too many dimensions.");
+                    bitsPerDimension = 64 / dimensionCount;
+                    if (dimensionCount == 1) {
+                        dimensionBitMask = -1ull;
+                    } else {
+                        dimensionBitMask = (1ull << bitsPerDimension) - 1;
+                    }
+                    
+                    if (dimensionCount * bitsPerDimension == 64ull) {
+                        relevantBitsMask = -1ull;
+                    } else {
+                        relevantBitsMask  = (1ull << (dimensionCount * bitsPerDimension)) - 1;
+                    }
+                }
+                
+                typename EpochManager::Epoch EpochManager::getBottomEpoch() const {
+                    return relevantBitsMask;
+                }
+                
+                typename EpochManager::Epoch EpochManager::getZeroEpoch() const {
+                    return 0;
+                }
+    
+                
+                uint64_t const& EpochManager::getDimensionCount() const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return dimensionCount;
+                }
+    
+                bool EpochManager::compareEpochClass(Epoch const& epoch1, Epoch const& epoch2) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    uint64_t mask = dimensionBitMask;
+                    for (uint64_t d = 0; d < dimensionCount; ++d) {
+                        if (((epoch1 & mask) == mask) != ((epoch2 & mask) == mask)) {
+                            assert(getEpochClass(epoch1) != getEpochClass(epoch2));
+                            return false;
+                        }
+                        mask = mask << bitsPerDimension;
+                    }
+                    assert(getEpochClass(epoch1) == getEpochClass(epoch2));
+                    return true;
+                }
+                
+                typename EpochManager::EpochClass EpochManager::getEpochClass(Epoch const& epoch) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    EpochClass result = 0;
+                    uint64_t dimensionMask = dimensionBitMask;
+                    uint64_t classMask = 1;
+                    for (uint64_t d  = 0; d < dimensionCount; ++d) {
+                        if ((epoch & dimensionMask) == dimensionMask) {
+                            result |= classMask;
+                        }
+                        classMask = classMask << 1;
+                        dimensionMask = dimensionMask << bitsPerDimension;
+                    }
+                    return result;
+                }
+    
+                typename EpochManager::Epoch EpochManager::getSuccessorEpoch(Epoch const& epoch, Epoch const& step) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    STORM_LOG_ASSERT(!hasBottomDimension(step), "The given step has at least one bottom dimension.");
+    
+                    // Start with dimension zero
+                    uint64_t mask = dimensionBitMask;
+                    uint64_t e_d = epoch & mask;
+                    uint64_t s_d = step & mask;
+                    uint64_t result = (e_d < s_d || e_d == mask) ? mask : e_d - s_d;
+                    
+                    // Consider the remaining dimensions
+                    for (uint64_t d = 1; d < dimensionCount; ++d) {
+                        mask = mask << bitsPerDimension;
+                        e_d = epoch & mask;
+                        s_d = step & mask;
+                        result |= ((e_d < s_d || e_d == mask) ? mask : e_d - s_d);
+                    }
+                    return result;
+                }
+                
+                std::vector<typename EpochManager::Epoch> EpochManager::getPredecessorEpochs(Epoch const& epoch, Epoch const& step) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    STORM_LOG_ASSERT(!hasBottomDimension(step), "The given step has at least one bottom dimension.");
+                    std::set<Epoch> resultAsSet;
+                    gatherPredecessorEpochs(resultAsSet, epoch, step);
+                    return std::vector<Epoch>(resultAsSet.begin(), resultAsSet.end());
+                }
+                
+                void EpochManager::gatherPredecessorEpochs(std::set<Epoch>& gatheredPredecessorEpochs, Epoch const& epoch, Epoch const& step) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    STORM_LOG_ASSERT(!hasBottomDimension(step), "The given step has at least one bottom dimension.");
+                    Epoch currStep = step;
+                    uint64_t d = 0;
+                    while (d < dimensionCount) {
+                        Epoch predecessor = epoch;
+                        for (uint64_t dPrime = 0; dPrime < dimensionCount; ++dPrime) {
+                            uint64_t step_dPrime = getDimensionOfEpoch(currStep, dPrime);
+                            if (isBottomDimension(predecessor, dPrime)) {
+                                if (step_dPrime != 0) {
+                                    setDimensionOfEpoch(predecessor, dPrime, step_dPrime - 1);
+                                }
+                            } else {
+                                setDimensionOfEpoch(predecessor, dPrime, getDimensionOfEpoch(predecessor, dPrime) + step_dPrime);
+                            }
+                        }
+                        assert(epoch == getSuccessorEpoch(predecessor, step));
+                        gatheredPredecessorEpochs.insert(predecessor);
+                        
+                        do {
+                            if (isBottomDimension(epoch, d)) {
+                                uint64_t step_d = getDimensionOfEpoch(currStep, d);
+                                if (step_d == 0) {
+                                    setDimensionOfEpoch(currStep, d, getDimensionOfEpoch(step, d));
+                                } else {
+                                    setDimensionOfEpoch(currStep, d, step_d - 1);
+                                    d = 0;
+                                    break;
+                                }
+                            }
+                            ++d;
+                        } while(d < dimensionCount);
+                    }
+                }
+                
+                bool EpochManager::isValidDimensionValue(uint64_t const& value) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return ((value & dimensionBitMask) == value) && value != dimensionBitMask;
+                }
+                
+                bool EpochManager::isZeroEpoch(Epoch const& epoch) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return (epoch & relevantBitsMask) == 0;
+                }
+                
+                bool EpochManager::isBottomEpoch(Epoch const& epoch) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return (epoch & relevantBitsMask) == relevantBitsMask;
+                }
+                
+                bool EpochManager::hasBottomDimension(Epoch const& epoch) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    uint64_t mask = dimensionBitMask;
+                    for (uint64_t d  = 0; d < dimensionCount; ++d) {
+                        if ((epoch | mask) == epoch) {
+                            return true;
+                        }
+                        mask = mask << bitsPerDimension;
+                    }
+                    return false;
+                }
+                
+                bool EpochManager::hasBottomDimensionEpochClass(EpochClass const& epochClass) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    uint64_t mask = (1 << dimensionCount) - 1;
+                    return (epochClass & mask) != 0;
+                }
+                
+                bool EpochManager::isPredecessorEpochClass(EpochClass const& epochClass1, EpochClass const& epochClass2) const {
+                    return (epochClass1 & epochClass2) == epochClass1;
+                }
+    
+                void EpochManager::setBottomDimension(Epoch& epoch, uint64_t const& dimension) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    epoch |= (dimensionBitMask << (dimension * bitsPerDimension));
+                }
+                
+                void EpochManager::setDimensionOfEpoch(Epoch& epoch, uint64_t const& dimension, uint64_t const& value) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    STORM_LOG_ASSERT(isValidDimensionValue(value), "The dimension value " << value << " is too high.");
+                    epoch &= ~(dimensionBitMask << (dimension * bitsPerDimension));
+                    epoch |= (value << (dimension * bitsPerDimension));
+                }
+                
+                bool EpochManager::isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return (epoch | (dimensionBitMask << (dimension * bitsPerDimension))) == epoch;
+                }
+                
+                bool EpochManager::isBottomDimensionEpochClass(EpochClass const& epochClass, uint64_t const& dimension) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return (epochClass | (1 << dimension)) == epochClass;
+                }
+                
+                uint64_t EpochManager::getDimensionOfEpoch(Epoch const& epoch, uint64_t const& dimension) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    return (epoch >> (dimension * bitsPerDimension)) & dimensionBitMask;
+                }
+                
+                std::string EpochManager::toString(Epoch const& epoch) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    std::string res = "<" + (isBottomDimension(epoch, 0) ? "_" : std::to_string(getDimensionOfEpoch(epoch, 0)));
+                    for (uint64_t d = 1; d < dimensionCount; ++d) {
+                        res += ", ";
+                        res += (isBottomDimension(epoch, d) ? "_" : std::to_string(getDimensionOfEpoch(epoch, d)));
+                    }
+                    res += ">";
+                    return res;
+                }
+                
+                bool EpochManager::epochClassZigZagOrder(Epoch const& epoch1, Epoch const& epoch2) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    
+                    // Return true iff epoch 1 has to be computed before epoch 2
+                    
+                    if (!compareEpochClass(epoch1, epoch2)) {
+                        return epochClassOrder(getEpochClass(epoch1), getEpochClass(epoch2));
+                    }
+                    
+                    // check whether the sum of dimensions is the same
+                    uint64_t e1Sum = 0;
+                    uint64_t e2Sum = 0;
+                    for (uint64_t dim = 0; dim < dimensionCount; ++dim) {
+                        if (!isBottomDimension(epoch1, dim)) {
+                            assert(!isBottomDimension(epoch2, dim));
+                            e1Sum += getDimensionOfEpoch(epoch1, dim);
+                            e2Sum += getDimensionOfEpoch(epoch2, dim);
+                        }
+                    }
+                    if (e1Sum < e2Sum) {
+                        return true;
+                    } else if (e1Sum > e2Sum) {
+                        return false;
+                    }
+                    
+                    // find the first dimension where the epochs do not match.
+                    // if the sum is even, we search from left to right, otherwise from right to left
+                    bool sumEven = (e1Sum % 2) == 0;
+                    if (sumEven) {
+                        for (uint64_t dim = 0; dim < dimensionCount; ++dim) {
+                            uint64_t e1Value = getDimensionOfEpoch(epoch1, dim);
+                            uint64_t e2Value = getDimensionOfEpoch(epoch2, dim);
+                            if (e1Value < e2Value) {
+                                return true;
+                            } else if (e1Value > e2Value) {
+                                return false;
+                            }
+                        }
+                    } else {
+                        uint64_t dim = dimensionCount;
+                        while (dim > 0) {
+                            --dim;
+                            uint64_t e1Value = getDimensionOfEpoch(epoch1, dim);
+                            uint64_t e2Value = getDimensionOfEpoch(epoch2, dim);
+                            if (e1Value < e2Value) {
+                                return true;
+                            } else if (e1Value > e2Value) {
+                                return false;
+                            }
+                        }
+                    }
+                    
+                    // reaching this point means that the epochs are equal
+                    assert(epoch1 == epoch2);
+                    return false;
+                }
+                
+                
+                bool EpochManager::epochClassOrder(EpochClass const& epochClass1, EpochClass const& epochClass2) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    
+                    if (epochClass1 == epochClass2) {
+                        return false;
+                    }
+                    
+                    // Return true iff epochClass 1 is a successor class of epochClass 2
+                    
+                    // Check whether the number of bottom dimensions is not equal
+                    int64_t count;
+                    EpochClass ec = epochClass1;
+                    for (count = 0; ec; ++count) {
+                        ec &= ec - 1;
+                    }
+                    ec = epochClass2;
+                    for (; ec; --count) {
+                        ec &= ec - 1;
+                    }
+                    if (count > 0) {
+                        return true;
+                    } else if (count < 0) {
+                        return false;
+                    }
+        
+                    return epochClass1 < epochClass2;
+                }
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
@ -0,0 +1,63 @@
+#pragma once
+
+#include <vector>
+#include <set>
+#include <string>
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+    
+                class EpochManager {
+                public:
+                    
+                    typedef uint64_t Epoch; // The number of reward steps that are "left" for each dimension
+                    typedef uint64_t EpochClass; // The number of reward steps that are "left" for each dimension
+                    
+                    EpochManager();
+                    EpochManager(uint64_t dimensionCount);
+                    
+                    uint64_t const& getDimensionCount() const;
+                    
+                    Epoch getBottomEpoch() const;
+                    Epoch getZeroEpoch() const;
+                    
+                    bool compareEpochClass(Epoch const& epoch1, Epoch const& epoch2) const;
+                    EpochClass getEpochClass(Epoch const& epoch) const;
+                    
+                    Epoch getSuccessorEpoch(Epoch const& epoch, Epoch const& step) const;
+                    std::vector<Epoch> getPredecessorEpochs(Epoch const& epoch, Epoch const& step) const;
+                    void gatherPredecessorEpochs(std::set<Epoch>& gatheredPredecessorEpochs, Epoch const& epoch, Epoch const& step) const;
+                    
+                    bool isZeroEpoch(Epoch const& epoch) const;
+                    bool isBottomEpoch(Epoch const& epoch) const;
+                    bool hasBottomDimension(Epoch const& epoch) const;
+                    bool hasBottomDimensionEpochClass(EpochClass const& epochClass) const;
+                    bool isPredecessorEpochClass(EpochClass const& epochClass1, EpochClass const& epochClass2) const;
+                    
+                    bool isValidDimensionValue(uint64_t const& value) const;
+                    
+                    void setBottomDimension(Epoch& epoch, uint64_t const& dimension) const;
+                    void setDimensionOfEpoch(Epoch& epoch, uint64_t const& dimension, uint64_t const& value) const; // assumes that the value is valid, i.e., small enough
+                    
+                    bool isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const;
+                    bool isBottomDimensionEpochClass(EpochClass const& epochClass, uint64_t const& dimension) const;
+                    uint64_t getDimensionOfEpoch(Epoch const& epoch, uint64_t const& dimension) const; // assumes that the dimension is not bottom
+                    
+                    std::string toString(Epoch const& epoch) const;
+                    
+                    bool epochClassZigZagOrder(Epoch const& epoch1, Epoch const& epoch2) const;
+                    bool epochClassOrder(EpochClass const& epochClass1, EpochClass const& epochClass2) const;
+                    
+                private:
+    
+                    uint64_t dimensionCount;
+                    uint64_t bitsPerDimension;
+                    uint64_t dimensionBitMask;
+                    uint64_t relevantBitsMask;
+                };
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.cpp
@ -0,0 +1,90 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h"
+
+#include "storm/utility/macros.h"
+
+#include "storm/exceptions/IllegalArgumentException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+    
+                MemoryStateManager::MemoryStateManager(uint64_t dimensionCount) : dimensionCount(dimensionCount), dimensionBitMask(1ull), relevantBitsMask((1ull << dimensionCount) - 1), stateCount(dimensionBitMask << dimensionCount), dimensionsWithoutMemoryMask(0), upperMemoryStateBound(stateCount) {
+                    STORM_LOG_THROW(dimensionCount > 0, storm::exceptions::IllegalArgumentException, "Invoked MemoryStateManager with zero dimension count.");
+                    STORM_LOG_THROW(dimensionCount <= 64, storm::exceptions::IllegalArgumentException, "Invoked MemoryStateManager with too many dimensions.");
+                }
+                
+                uint64_t const& MemoryStateManager::getDimensionCount() const {
+                    return dimensionCount;
+                }
+                
+                uint64_t const& MemoryStateManager::getMemoryStateCount() const {
+                    return stateCount;
+                }
+                
+                MemoryStateManager::MemoryState const& MemoryStateManager::getUpperMemoryStateBound() const {
+                    return upperMemoryStateBound;
+                }
+                
+                void MemoryStateManager::setDimensionWithoutMemory(uint64_t dimension) {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked MemoryStateManager with zero dimension count.");
+                    STORM_LOG_ASSERT(dimension < dimensionCount, "Tried to set a dimension that is larger then the number of considered dimensions");
+                    if (((dimensionBitMask << dimension) & dimensionsWithoutMemoryMask) == 0) {
+                        stateCount = (stateCount >> 1);
+                    }
+                    dimensionsWithoutMemoryMask |= (dimensionBitMask << dimension);
+                }
+    
+                MemoryStateManager::MemoryState MemoryStateManager::getInitialMemoryState() const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked MemoryStateManager with zero dimension count.");
+                    return relevantBitsMask;
+                }
+                
+                bool MemoryStateManager::isRelevantDimension(MemoryState const& state, uint64_t dimension) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked MemoryStateManager with zero dimension count.");
+                    STORM_LOG_ASSERT((state & dimensionsWithoutMemoryMask) == dimensionsWithoutMemoryMask, "Invalid memory state found.");
+                    return (state & (dimensionBitMask << dimension)) != 0;
+                }
+                
+                void MemoryStateManager::setRelevantDimension(MemoryState& state, uint64_t dimension, bool value) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked MemoryStateManager with zero dimension count.");
+                    STORM_LOG_ASSERT(dimension < dimensionCount, "Tried to set a dimension that is larger then the number of considered dimensions");
+                    STORM_LOG_ASSERT(((dimensionBitMask << dimension) & dimensionsWithoutMemoryMask) == 0, "Tried to change a memory state for a dimension but the dimension is assumed to have no memory.");
+                    if (value) {
+                        state |= (dimensionBitMask << dimension);
+                    } else {
+                        state &= ~(dimensionBitMask << dimension);
+                    }
+                }
+                
+                void MemoryStateManager::setRelevantDimensions(MemoryState& state, storm::storage::BitVector const& dimensions, bool value) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked MemoryStateManager with zero dimension count.");
+                    STORM_LOG_ASSERT(dimensions.size() == dimensionCount, "Invalid size of given bitset.");
+                    if (value) {
+                        for (auto const& d : dimensions) {
+                            STORM_LOG_ASSERT(((dimensionBitMask << d) & dimensionsWithoutMemoryMask) == 0, "Tried to set a dimension to 'relevant'-memory state but the dimension is assumed to have no memory.");
+                            state |= (dimensionBitMask << d);
+                        }
+                    } else {
+                        for (auto const& d : dimensions) {
+                            STORM_LOG_ASSERT(((dimensionBitMask << d) & dimensionsWithoutMemoryMask) == 0, "Tried to set a dimension to 'unrelevant'-memory state but the dimension is assumed to have no memory.");
+                            state &= ~(dimensionBitMask << d);
+                        }
+                    }
+                }
+                
+                std::string MemoryStateManager::toString(MemoryState const& state) const {
+                    STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
+                    std::string res = "[";
+                    res += (isRelevantDimension(state, 0) ? "1" : "0");
+                    for (uint64_t d = 1; d < dimensionCount; ++d) {
+                        res += ", ";
+                        res += (isRelevantDimension(state, d) ? "1" : "0");
+                    }
+                    res += "]";
+                    return res;
+                }
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h
@ -0,0 +1,46 @@
+#pragma once
+
+#include <vector>
+#include <set>
+#include <string>
+
+#include "storm/storage/BitVector.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+            
+                class MemoryStateManager {
+                public:
+                    
+                    typedef uint64_t MemoryState;
+                    
+                    MemoryStateManager(uint64_t dimensionCount);
+                    
+                    void setDimensionWithoutMemory(uint64_t dimension);
+                    
+                    uint64_t const& getDimensionCount() const;
+                    uint64_t const& getMemoryStateCount() const;
+                    MemoryState const& getUpperMemoryStateBound() const; // is larger then every valid memory state m, i.e., m < getUpperMemoryStateBound() holds for all m
+                    
+                    MemoryState getInitialMemoryState() const;
+                    bool isRelevantDimension(MemoryState const& state, uint64_t dimension) const;
+                    void setRelevantDimension(MemoryState& state, uint64_t dimension, bool value = true) const;
+                    void setRelevantDimensions(MemoryState& state, storm::storage::BitVector const& dimensions, bool value = true) const;
+                    
+                    std::string toString(MemoryState const& state) const;
+                    
+                private:
+                    uint64_t const dimensionCount;
+                    uint64_t const dimensionBitMask;
+                    uint64_t const relevantBitsMask;
+                    uint64_t stateCount;
+                    uint64_t dimensionsWithoutMemoryMask;
+                    MemoryState const upperMemoryStateBound;
+    
+                };
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
@ -0,0 +1,640 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
+
+#include <string>
+#include <set>
+#include <functional>
+
+#include "storm/utility/macros.h"
+#include "storm/logic/Formulas.h"
+
+#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/transformer/EndComponentEliminator.h"
+
+#include "storm/exceptions/UnexpectedException.h"
+#include "storm/exceptions/IllegalArgumentException.h"
+#include "storm/exceptions/NotSupportedException.h"
+#include "storm/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+    
+                template<typename ValueType, bool SingleObjectiveMode>
+                MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) : model(model), objectives(objectives) {
+                    initialize();
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> objectiveFormula) : model(model) {
+                    
+                    STORM_LOG_THROW(objectiveFormula->hasOptimalityType(), storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+                    if (objectiveFormula->isProbabilityOperatorFormula()) {
+                        if (objectiveFormula->getSubformula().isMultiObjectiveFormula()) {
+                            for (auto const& subFormula : objectiveFormula->getSubformula().asMultiObjectiveFormula().getSubformulas()) {
+                                STORM_LOG_THROW(subFormula->isBoundedUntilFormula(), storm::exceptions::InvalidPropertyException, "Formula " << objectiveFormula << " is not supported. Invalid subformula " << *subFormula << ".");
+                            }
+                        } else {
+                            STORM_LOG_THROW(objectiveFormula->getSubformula().isBoundedUntilFormula(), storm::exceptions::InvalidPropertyException, "Formula " << objectiveFormula << " is not supported. Invalid subformula " << objectiveFormula->getSubformula() << ".");
+                        }
+                    } else {
+                        STORM_LOG_THROW(objectiveFormula->isRewardOperatorFormula() && objectiveFormula->getSubformula().isCumulativeRewardFormula(), storm::exceptions::InvalidPropertyException, "Formula " << objectiveFormula << " is not supported.");
+                        
+                    }
+                    
+                    // Build an objective from the formula.
+                    storm::modelchecker::multiobjective::Objective<ValueType> objective;
+                    objective.formula = objectiveFormula;
+                    objective.originalFormula = objective.formula;
+                    objective.considersComplementaryEvent = false;
+                    objectives.push_back(std::move(objective));
+                    
+                    initialize();
+                }
+        
+        
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initialize() {
+                    
+                    maxSolutionsStored = 0;
+                    
+                    STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives.");
+                    std::vector<Epoch> epochSteps;
+                    initializeObjectives(epochSteps);
+                    initializeMemoryProduct(epochSteps);
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initializeObjectives(std::vector<Epoch>& epochSteps) {
+                    std::vector<std::vector<uint64_t>> dimensionWiseEpochSteps;
+                    // collect the time-bounded subobjectives
+                    for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                        auto const& formula = *this->objectives[objIndex].formula;
+                        if (formula.isProbabilityOperatorFormula()) {
+                            STORM_LOG_THROW(formula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Unexpected type of subformula for formula " << formula);
+                            auto const& subformula = formula.getSubformula().asBoundedUntilFormula();
+                            for (uint64_t dim = 0; dim < subformula.getDimension(); ++dim) {
+                                Dimension<ValueType> dimension;
+                                dimension.formula = subformula.restrictToDimension(dim);
+                                dimension.objectiveIndex = objIndex;
+                                std::string memLabel = "dim" + std::to_string(dimensions.size()) + "_maybe";
+                                while (model.getStateLabeling().containsLabel(memLabel)) {
+                                    memLabel = "_" + memLabel;
+                                }
+                                dimension.memoryLabel = memLabel;
+                                dimension.isUpperBounded = subformula.hasUpperBound(dim);
+                                // for simplicity we do not allow intervals or unbounded formulas.
+                                STORM_LOG_THROW(subformula.hasLowerBound(dim) != dimension.isUpperBounded, storm::exceptions::NotSupportedException, "Bounded until formulas are only supported by this method if they consider either an upper bound or a lower bound. Got " << subformula << " instead.");
+                                // lower bounded until formulas with non-trivial left hand side are excluded as this would require some additional effort (in particular the ProductModel::transformMemoryState method).
+                                STORM_LOG_THROW(dimension.isUpperBounded || subformula.getLeftSubformula(dim).isTrueFormula(), storm::exceptions::NotSupportedException, "Lower bounded until formulas are only supported by this method if the left subformula is 'true'. Got " << subformula << " instead.");
+                                if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) {
+                                    dimensionWiseEpochSteps.push_back(std::vector<uint64_t>(model.getNumberOfChoices(), 1));
+                                    dimension.scalingFactor = storm::utility::one<ValueType>();
+                                } else {
+                                    STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound.");
+                                    std::string const& rewardName = subformula.getTimeBoundReference(dim).getRewardName();
+                                    STORM_LOG_THROW(this->model.hasRewardModel(rewardName), storm::exceptions::IllegalArgumentException, "No reward model with name '" << rewardName << "' found.");
+                                    auto const& rewardModel = this->model.getRewardModel(rewardName);
+                                    STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Transition rewards are currently not supported as reward bounds.");
+                                    std::vector<ValueType> actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix());
+                                    auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector<ValueType, uint64_t>(actionRewards);
+                                    dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first));
+                                    dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second);
+                                }
+                                dimensions.emplace_back(std::move(dimension));
+                            }
+                        } else if (formula.isRewardOperatorFormula() && formula.getSubformula().isCumulativeRewardFormula()) {
+                            auto const& subformula = formula.getSubformula().asCumulativeRewardFormula();
+                            Dimension<ValueType> dimension;
+                            dimension.formula = formula.getSubformula().asSharedPointer();
+                            dimension.objectiveIndex = objIndex;
+                            dimension.isUpperBounded = true;
+                            if (subformula.getTimeBoundReference().isTimeBound() || subformula.getTimeBoundReference().isStepBound()) {
+                                dimensionWiseEpochSteps.push_back(std::vector<uint64_t>(model.getNumberOfChoices(), 1));
+                                dimension.scalingFactor = storm::utility::one<ValueType>();
+                            } else {
+                                STORM_LOG_ASSERT(subformula.getTimeBoundReference().isRewardBound(), "Unexpected type of time bound.");
+                                std::string const& rewardName = subformula.getTimeBoundReference().getRewardName();
+                                STORM_LOG_THROW(this->model.hasRewardModel(rewardName), storm::exceptions::IllegalArgumentException, "No reward model with name '" << rewardName << "' found.");
+                                auto const& rewardModel = this->model.getRewardModel(rewardName);
+                                STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Transition rewards are currently not supported as reward bounds.");
+                                std::vector<ValueType> actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix());
+                                auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector<ValueType, uint64_t>(actionRewards);
+                                dimensionWiseEpochSteps.push_back(std::move(discretizedRewardsAndFactor.first));
+                                dimension.scalingFactor = std::move(discretizedRewardsAndFactor.second);
+                            }
+                            dimensions.emplace_back(std::move(dimension));
+                        }
+                    }
+                    
+                    // Compute a mapping for each objective to the set of dimensions it considers
+                    // Also store for each dimension dim, which dimensions should be unsatisfiable whenever the bound of dim is violated
+                    uint64_t dim = 0;
+                    for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                        storm::storage::BitVector objDimensions(dimensions.size(), false);
+                        if (objectives[objIndex].formula->isProbabilityOperatorFormula() && objectives[objIndex].formula->getSubformula().isBoundedUntilFormula()) {
+                            auto const& boundedUntilFormula = objectives[objIndex].formula->getSubformula().asBoundedUntilFormula();
+                            for (uint64_t currDim = dim; currDim < dim + boundedUntilFormula.getDimension(); ++currDim ) {
+                                objDimensions.set(currDim);
+                            }
+                            for (uint64_t currDim = dim; currDim < dim + boundedUntilFormula.getDimension(); ++currDim ) {
+                                if (!boundedUntilFormula.hasMultiDimensionalSubformulas() || dimensions[currDim].isUpperBounded) {
+                                    dimensions[currDim].dependentDimensions = objDimensions;
+                                } else {
+                                    dimensions[currDim].dependentDimensions = storm::storage::BitVector(dimensions.size(), false);
+                                    dimensions[currDim].dependentDimensions.set(currDim, true);
+                                }
+                                //  std::cout << "dimension " << currDim << " has depDims: " << dimensions[currDim].dependentDimensions << std::endl;
+                            }
+                            dim += boundedUntilFormula.getDimension();
+                        } else if (objectives[objIndex].formula->isRewardOperatorFormula() && objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) {
+                            objDimensions.set(dim, true);
+                            dimensions[dim].dependentDimensions = objDimensions;
+                            ++dim;
+                        }
+                        
+                        objectiveDimensions.push_back(std::move(objDimensions));
+                    }
+                    assert(dim == dimensions.size());
+                    
+                    
+                    // Initialize the epoch manager
+                    epochManager = EpochManager(dimensions.size());
+                    
+                    // Convert the epoch steps to a choice-wise representation
+                    epochSteps.reserve(model.getNumberOfChoices());
+                    for (uint64_t choice = 0; choice < model.getNumberOfChoices(); ++choice) {
+                        Epoch step;
+                        uint64_t dim = 0;
+                        for (auto const& dimensionSteps : dimensionWiseEpochSteps) {
+                            epochManager.setDimensionOfEpoch(step, dim, dimensionSteps[choice]);
+                            ++dim;
+                        }
+                        epochSteps.push_back(step);
+                    }
+                    
+                    // collect which epoch steps are possible
+                    possibleEpochSteps.clear();
+                    for (auto const& step : epochSteps) {
+                        possibleEpochSteps.insert(step);
+                    }
+                    
+                    // Set the maximal values we need to consider for each dimension
+                    computeMaxDimensionValues();
+                    
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initializeMemoryProduct(std::vector<Epoch> const& epochSteps) {
+                    productModel = std::make_unique<ProductModel<ValueType>>(model, objectives, dimensions, objectiveDimensions, epochManager, epochSteps);
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::computeMaxDimensionValues() {
+                    for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                        storm::expressions::Expression bound;
+                        bool isStrict = false;
+                        storm::logic::Formula const& dimFormula = *dimensions[dim].formula;
+                        if (dimFormula.isBoundedUntilFormula()) {
+                            assert(!dimFormula.asBoundedUntilFormula().isMultiDimensional());
+                            if (dimFormula.asBoundedUntilFormula().hasUpperBound()) {
+                                STORM_LOG_ASSERT(!dimFormula.asBoundedUntilFormula().hasLowerBound(), "Bounded until formulas with interval bounds are not supported.");
+                                bound = dimFormula.asBoundedUntilFormula().getUpperBound();
+                                isStrict = dimFormula.asBoundedUntilFormula().isUpperBoundStrict();
+                            } else {
+                                STORM_LOG_ASSERT(dimFormula.asBoundedUntilFormula().hasLowerBound(), "Bounded until formulas without any bounds are not supported.");
+                                bound = dimFormula.asBoundedUntilFormula().getLowerBound();
+                                isStrict = dimFormula.asBoundedUntilFormula().isLowerBoundStrict();
+                            }
+                        } else if (dimFormula.isCumulativeRewardFormula()) {
+                            bound = dimFormula.asCumulativeRewardFormula().getBound();
+                            isStrict = dimFormula.asCumulativeRewardFormula().isBoundStrict();
+                        }
+                        STORM_LOG_THROW(!bound.containsVariables(), storm::exceptions::NotSupportedException, "The bound " << bound << " contains undefined constants.");
+                        ValueType discretizedBound = storm::utility::convertNumber<ValueType>(bound.evaluateAsRational());
+                        STORM_LOG_THROW(dimensions[dim].isUpperBounded || isStrict || !storm::utility::isZero(discretizedBound), storm::exceptions::NotSupportedException, "Lower bounds need to be either strict or greater than zero.");
+                        discretizedBound /= dimensions[dim].scalingFactor;
+                        if (storm::utility::isInteger(discretizedBound)) {
+                            if (isStrict == dimensions[dim].isUpperBounded) {
+                                discretizedBound -= storm::utility::one<ValueType>();
+                            }
+                        } else {
+                            discretizedBound = storm::utility::floor(discretizedBound);
+                        }
+                        uint64_t dimensionValue = storm::utility::convertNumber<uint64_t>(discretizedBound);
+                        STORM_LOG_THROW(epochManager.isValidDimensionValue(dimensionValue), storm::exceptions::NotSupportedException, "The bound " << bound << " is too high for the considered number of dimensions.");
+                        dimensions[dim].maxValue = dimensionValue;
+                    }
+                }
+        
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStartEpoch() {
+                    Epoch startEpoch = epochManager.getZeroEpoch();
+                    for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                        STORM_LOG_ASSERT(dimensions[dim].maxValue,  "No max-value for dimension " << dim << " was given.");
+                        epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get());
+                    }
+                    STORM_LOG_TRACE("Start epoch is " << epochManager.toString(startEpoch));
+                    return startEpoch;
+                }
+        
+                template<typename ValueType, bool SingleObjectiveMode>
+                std::vector<typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochComputationOrder(Epoch const& startEpoch) {
+                    // Perform a DFS to find all the reachable epochs
+                    std::vector<Epoch> dfsStack;
+                    std::set<Epoch, std::function<bool(Epoch const&, Epoch const&)>> collectedEpochs(std::bind(&EpochManager::epochClassZigZagOrder, &epochManager, std::placeholders::_1, std::placeholders::_2));
+                    
+                    collectedEpochs.insert(startEpoch);
+                    dfsStack.push_back(startEpoch);
+                    while (!dfsStack.empty()) {
+                        Epoch currentEpoch = dfsStack.back();
+                        dfsStack.pop_back();
+                        for (auto const& step : possibleEpochSteps) {
+                            Epoch successorEpoch = epochManager.getSuccessorEpoch(currentEpoch, step);
+                            /*
+                            for (auto const& e : collectedEpochs) {
+                                std::cout << "Comparing " << epochManager.toString(e) << " and " << epochManager.toString(successorEpoch) << std::endl;
+                                if (epochManager.epochClassZigZagOrder(e, successorEpoch)) {
+                                    std::cout << "    " << epochManager.toString(e) << " < " << epochManager.toString(successorEpoch) << std::endl;
+                                }
+                                if (epochManager.epochClassZigZagOrder(successorEpoch, e)) {
+                                    std::cout << "    " << epochManager.toString(e) << " > " << epochManager.toString(successorEpoch) << std::endl;
+                                }
+                            }
+                             */
+                            if (collectedEpochs.insert(successorEpoch).second) {
+                                dfsStack.push_back(std::move(successorEpoch));
+                            }
+                        }
+                    }
+                    /*
+                    std::cout << "Resulting order: ";
+                    for (auto const& e : collectedEpochs) {
+                        std::cout << epochManager.toString(e) << ", ";
+                    }
+                    std::cout << std::endl;
+                    */
+                    return std::vector<Epoch>(collectedEpochs.begin(), collectedEpochs.end());
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::EpochModel& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpoch(Epoch const& epoch) {
+                    STORM_LOG_DEBUG("Setting model for epoch " << epochManager.toString(epoch));
+                    
+                    // Check if we need to update the current epoch class
+                    if (!currentEpoch || !epochManager.compareEpochClass(epoch, currentEpoch.get())) {
+                        setCurrentEpochClass(epoch);
+                        epochModel.epochMatrixChanged = true;
+                    } else {
+                        epochModel.epochMatrixChanged = false;
+                    }
+                    
+                    bool containsLowerBoundedObjective = false;
+                    for (auto const& dimension : dimensions) {
+                        if (!dimension.isUpperBounded) {
+                            containsLowerBoundedObjective = true;
+                            break;
+                        }
+                    }
+                    std::map<Epoch, EpochSolution const*> subSolutions;
+                    for (auto const& step : possibleEpochSteps) {
+                        Epoch successorEpoch = epochManager.getSuccessorEpoch(epoch, step);
+                        if (successorEpoch != epoch) {
+                            auto successorSolIt = epochSolutions.find(successorEpoch);
+                            STORM_LOG_ASSERT(successorSolIt != epochSolutions.end(), "Solution for successor epoch does not exist (anymore).");
+                            subSolutions.emplace(successorEpoch, &successorSolIt->second);
+                        }
+                    }
+                    epochModel.stepSolutions.resize(epochModel.stepChoices.getNumberOfSetBits());
+                    auto stepSolIt = epochModel.stepSolutions.begin();
+                    for (auto const& reducedChoice : epochModel.stepChoices) {
+                        uint64_t productChoice = epochModelToProductChoiceMap[reducedChoice];
+                        uint64_t productState = productModel->getProductStateFromChoice(productChoice);
+                        auto const& memoryState = productModel->getMemoryState(productState);
+                        Epoch successorEpoch = epochManager.getSuccessorEpoch(epoch, productModel->getSteps()[productChoice]);
+                        EpochClass successorEpochClass = epochManager.getEpochClass(successorEpoch);
+                        // Find out whether objective reward is earned for the current choice
+                        // Objective reward is not earned if
+                        // a) there is an upper bounded subObjective that is __still_relevant__ but the corresponding reward bound is passed after taking the choice
+                        // b) there is a lower bounded subObjective and the corresponding reward bound is not passed yet.
+                        for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                            bool rewardEarned = !storm::utility::isZero(epochModel.objectiveRewards[objIndex][reducedChoice]);
+                            if (rewardEarned) {
+                                for (auto const& dim : objectiveDimensions[objIndex]) {
+                                    if (dimensions[dim].isUpperBounded == epochManager.isBottomDimension(successorEpoch, dim) && productModel->getMemoryStateManager().isRelevantDimension(memoryState, dim)) {
+                                        rewardEarned = false;
+                                        break;
+                                    }
+                                }
+                            }
+                            epochModel.objectiveRewardFilter[objIndex].set(reducedChoice, rewardEarned);
+                        }
+                        // compute the solution for the stepChoices
+                        // For optimization purposes, we distinguish the case where the memory state does not have to be transformed
+                        EpochSolution const& successorEpochSolution = getEpochSolution(subSolutions, successorEpoch);
+                        SolutionType choiceSolution;
+                        bool firstSuccessor = true;
+                        if (!containsLowerBoundedObjective && epochManager.compareEpochClass(epoch, successorEpoch)) {
+                            for (auto const& successor : productModel->getProduct().getTransitionMatrix().getRow(productChoice)) {
+                                if (firstSuccessor) {
+                                    choiceSolution = getScaledSolution(getStateSolution(successorEpochSolution, successor.getColumn()), successor.getValue());
+                                    firstSuccessor = false;
+                                } else {
+                                    addScaledSolution(choiceSolution, getStateSolution(successorEpochSolution, successor.getColumn()), successor.getValue());
+                                }
+                            }
+                        } else {
+                            for (auto const& successor : productModel->getProduct().getTransitionMatrix().getRow(productChoice)) {
+                                uint64_t successorProductState = productModel->transformProductState(successor.getColumn(), successorEpochClass, memoryState);
+                                SolutionType const& successorSolution = getStateSolution(successorEpochSolution, successorProductState);
+                                if (firstSuccessor) {
+                                    choiceSolution = getScaledSolution(successorSolution, successor.getValue());
+                                    firstSuccessor = false;
+                                } else {
+                                    addScaledSolution(choiceSolution, successorSolution, successor.getValue());
+                                }
+                            }
+                        }
+                        assert(!firstSuccessor);
+                        *stepSolIt = std::move(choiceSolution);
+                        ++stepSolIt;
+                    }
+                    
+                    assert(epochModel.objectiveRewards.size() == objectives.size());
+                    assert(epochModel.objectiveRewardFilter.size() == objectives.size());
+                    assert(epochModel.epochMatrix.getRowCount() == epochModel.stepChoices.size());
+                    assert(epochModel.stepChoices.size() == epochModel.objectiveRewards.front().size());
+                    assert(epochModel.objectiveRewards.front().size() == epochModel.objectiveRewards.back().size());
+                    assert(epochModel.objectiveRewards.front().size() == epochModel.objectiveRewardFilter.front().size());
+                    assert(epochModel.objectiveRewards.back().size() == epochModel.objectiveRewardFilter.back().size());
+                    assert(epochModel.stepChoices.getNumberOfSetBits() == epochModel.stepSolutions.size());
+                    
+                    currentEpoch = epoch;
+                    /*
+                    std::cout << "Epoch model for epoch " << storm::utility::vector::toString(epoch) << std::endl;
+                    std::cout << "Matrix: " << std::endl << epochModel.epochMatrix << std::endl;
+                    std::cout << "ObjectiveRewards: " << storm::utility::vector::toString(epochModel.objectiveRewards[0]) << std::endl;
+                    std::cout << "steps: " << epochModel.stepChoices << std::endl;
+                    std::cout << "step solutions: ";
+                    for (int i = 0; i < epochModel.stepSolutions.size(); ++i) {
+                        std::cout << "   " << epochModel.stepSolutions[i].weightedValue;
+                    }
+                    std::cout << std::endl;
+                    */
+                    return epochModel;
+                    
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpochClass(Epoch const& epoch) {
+                    EpochClass epochClass = epochManager.getEpochClass(epoch);
+                    // std::cout << "Setting epoch class for epoch " << epochManager.toString(epoch) << std::endl;
+                    auto productObjectiveRewards = productModel->computeObjectiveRewards(epochClass, objectives);
+                    
+                    storm::storage::BitVector stepChoices(productModel->getProduct().getNumberOfChoices(), false);
+                    uint64_t choice = 0;
+                    for (auto const& step : productModel->getSteps()) {
+                        if (!epochManager.isZeroEpoch(step) && epochManager.getSuccessorEpoch(epoch, step) != epoch) {
+                            stepChoices.set(choice, true);
+                        }
+                        ++choice;
+                    }
+                    epochModel.epochMatrix = productModel->getProduct().getTransitionMatrix().filterEntries(~stepChoices);
+                    // redirect transitions for the case where the lower reward bounds are not met yet
+                    storm::storage::BitVector violatedLowerBoundedDimensions(dimensions.size(), false);
+                    for (uint64_t dim = 0; dim < dimensions.size(); ++dim) {
+                        if (!dimensions[dim].isUpperBounded && !epochManager.isBottomDimensionEpochClass(epochClass, dim)) {
+                            violatedLowerBoundedDimensions.set(dim);
+                        }
+                    }
+                    if (!violatedLowerBoundedDimensions.empty()) {
+                        for (uint64_t state = 0; state < epochModel.epochMatrix.getRowGroupCount(); ++state) {
+                            auto const& memoryState = productModel->getMemoryState(state);
+                            for (auto& entry : epochModel.epochMatrix.getRowGroup(state)) {
+                                entry.setColumn(productModel->transformProductState(entry.getColumn(), epochClass, memoryState));
+                            }
+                        }
+                    }
+                    
+                    storm::storage::BitVector zeroObjRewardChoices(productModel->getProduct().getNumberOfChoices(), true);
+                    for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                        if (violatedLowerBoundedDimensions.isDisjointFrom(objectiveDimensions[objIndex])) {
+                            zeroObjRewardChoices &= storm::utility::vector::filterZero(productObjectiveRewards[objIndex]);
+                        }
+                    }
+                    storm::storage::BitVector allProductStates(productModel->getProduct().getNumberOfStates(), true);
+                    
+                    // Get the relevant states for this epoch.
+                    storm::storage::BitVector productInStates = productModel->getInStates(epochClass);
+                    // The epoch model only needs to consider the states that are reachable from a relevant state
+                    storm::storage::BitVector consideredStates = storm::utility::graph::getReachableStates(epochModel.epochMatrix, productInStates, allProductStates, ~allProductStates);
+                    // std::cout << "numInStates = " << productInStates.getNumberOfSetBits() << std::endl;
+                    // std::cout << "numConsideredStates = " << consideredStates.getNumberOfSetBits() << std::endl;
+                    
+                    // We assume that there is no end component in which objective reward is earned
+                    STORM_LOG_ASSERT(!storm::utility::graph::checkIfECWithChoiceExists(epochModel.epochMatrix, epochModel.epochMatrix.transpose(true), allProductStates, ~zeroObjRewardChoices & ~stepChoices), "There is a scheduler that yields infinite reward for one objective. This case should be excluded");
+                    auto ecElimResult = storm::transformer::EndComponentEliminator<ValueType>::transform(epochModel.epochMatrix, consideredStates, zeroObjRewardChoices & ~stepChoices, consideredStates);
+                    epochModel.epochMatrix = std::move(ecElimResult.matrix);
+                    epochModelToProductChoiceMap = std::move(ecElimResult.newToOldRowMapping);
+                    
+                    epochModel.stepChoices = storm::storage::BitVector(epochModel.epochMatrix.getRowCount(), false);
+                    for (uint64_t choice = 0; choice < epochModel.epochMatrix.getRowCount(); ++choice) {
+                        if (stepChoices.get(epochModelToProductChoiceMap[choice])) {
+                            epochModel.stepChoices.set(choice, true);
+                        }
+                    }
+                    
+                    epochModel.objectiveRewards.clear();
+                    for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                        std::vector<ValueType> const& productObjRew = productObjectiveRewards[objIndex];
+                        std::vector<ValueType> reducedModelObjRewards;
+                        reducedModelObjRewards.reserve(epochModel.epochMatrix.getRowCount());
+                        for (auto const& productChoice : epochModelToProductChoiceMap) {
+                            reducedModelObjRewards.push_back(productObjRew[productChoice]);
+                        }
+                        // Check if the objective is violated in the current epoch
+                        if (!violatedLowerBoundedDimensions.isDisjointFrom(objectiveDimensions[objIndex])) {
+                            storm::utility::vector::setVectorValues(reducedModelObjRewards, ~epochModel.stepChoices, storm::utility::zero<ValueType>());
+                        }
+                        epochModel.objectiveRewards.push_back(std::move(reducedModelObjRewards));
+                    }
+                    
+                    epochModel.epochInStates = storm::storage::BitVector(epochModel.epochMatrix.getRowGroupCount(), false);
+                    for (auto const& productState : productInStates) {
+                        STORM_LOG_ASSERT(ecElimResult.oldToNewStateMapping[productState] < epochModel.epochMatrix.getRowGroupCount(), "Selected product state does not exist in the epoch model.");
+                        epochModel.epochInStates.set(ecElimResult.oldToNewStateMapping[productState], true);
+                    }
+                    
+                    std::vector<uint64_t> toEpochModelInStatesMap(productModel->getProduct().getNumberOfStates(), std::numeric_limits<uint64_t>::max());
+                    std::vector<uint64_t> epochModelStateToInStateMap = epochModel.epochInStates.getNumberOfSetBitsBeforeIndices();
+                    for (auto const& productState : productInStates) {
+                        toEpochModelInStatesMap[productState] = epochModelStateToInStateMap[ecElimResult.oldToNewStateMapping[productState]];
+                    }
+                    productStateToEpochModelInStateMap = std::make_shared<std::vector<uint64_t> const>(std::move(toEpochModelInStatesMap));
+                    
+                    epochModel.objectiveRewardFilter.clear();
+                    for (auto const& objRewards : epochModel.objectiveRewards) {
+                        epochModel.objectiveRewardFilter.push_back(storm::utility::vector::filterZero(objRewards));
+                        epochModel.objectiveRewardFilter.back().complement();
+                    }
+    
+                    epochModelSizes.push_back(epochModel.epochMatrix.getRowGroupCount());
+                }
+                
+     
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<SO, int>::type>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const {
+                    return solution * scalingFactor;
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<!SO, int>::type>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const {
+                    SolutionType res;
+                    res.reserve(solution.size());
+                    for (auto const& sol : solution) {
+                        res.push_back(sol * scalingFactor);
+                    }
+                    return res;
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<SO, int>::type>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const {
+                    solution += solutionToAdd * scalingFactor;
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<!SO, int>::type>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const {
+                    storm::utility::vector::addScaledVector(solution,  solutionToAdd, scalingFactor);
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<SO, int>::type>
+                std::string MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::solutionToString(SolutionType const& solution) const {
+                    std::stringstream stringstream;
+                    stringstream << solution;
+                    return stringstream.str();
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                template<bool SO, typename std::enable_if<!SO, int>::type>
+                std::string MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::solutionToString(SolutionType const& solution) const {
+                    std::stringstream stringstream;
+                    stringstream << "(";
+                    bool first = true;
+                    for (auto const& s : solution) {
+                        if (first) {
+                            first = false;
+                        } else {
+                            stringstream << ", ";
+                        }
+                        stringstream << s;
+                    }
+                    stringstream << ")";
+                    return stringstream.str();
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                ValueType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) {
+                    // if (storm::settings::getModule<storm::settings::modules::GeneralSettings>().isSoundSet()) {
+                    uint64_t sumOfDimensions = 0;
+                    for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                        sumOfDimensions += epochManager.getDimensionOfEpoch(startEpoch, dim) + 1;
+                    }
+                    return precision / storm::utility::convertNumber<ValueType>(sumOfDimensions);
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions) {
+                    STORM_LOG_ASSERT(currentEpoch, "Tried to set a solution for the current epoch, but no epoch was specified before.");
+                    STORM_LOG_ASSERT(inStateSolutions.size() == epochModel.epochInStates.getNumberOfSetBits(), "Invalid number of solutions.");
+    
+                    std::set<Epoch> predecessorEpochs, successorEpochs;
+                    for (auto const& step : possibleEpochSteps) {
+                        epochManager.gatherPredecessorEpochs(predecessorEpochs, currentEpoch.get(), step);
+                        successorEpochs.insert(epochManager.getSuccessorEpoch(currentEpoch.get(), step));
+                    }
+                    predecessorEpochs.erase(currentEpoch.get());
+                    successorEpochs.erase(currentEpoch.get());
+                    STORM_LOG_ASSERT(!predecessorEpochs.empty(), "There are no predecessors for the epoch " << epochManager.toString(currentEpoch.get()));
+                    
+                    // clean up solutions that are not needed anymore
+                    for (auto const& successorEpoch : successorEpochs) {
+                        auto successorEpochSolutionIt = epochSolutions.find(successorEpoch);
+                        STORM_LOG_ASSERT(successorEpochSolutionIt != epochSolutions.end(), "Solution for successor epoch does not exist (anymore).");
+                        --successorEpochSolutionIt->second.count;
+                        if (successorEpochSolutionIt->second.count == 0) {
+                            epochSolutions.erase(successorEpochSolutionIt);
+                        }
+                    }
+                    
+                    // add the new solution
+                    EpochSolution solution;
+                    solution.count = predecessorEpochs.size();
+                    solution.productStateToSolutionVectorMap = productStateToEpochModelInStateMap;
+                    solution.solutions = std::move(inStateSolutions);
+                    epochSolutions[currentEpoch.get()] = std::move(solution);
+                    
+                    maxSolutionsStored = std::max((uint64_t) epochSolutions.size(), maxSolutionsStored);
+                    
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStateSolution(Epoch const& epoch, uint64_t const& productState) {
+                    auto epochSolutionIt = epochSolutions.find(epoch);
+                    STORM_LOG_ASSERT(epochSolutionIt != epochSolutions.end(), "Requested unexisting solution for epoch " << epochManager.toString(epoch) << ".");
+                    auto const& epochSolution = epochSolutionIt->second;
+                    STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution for epoch " << epochManager.toString(epoch) << " at an unexisting product state.");
+                    STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch " << epochManager.toString(epoch) << " at a state for which no solution was stored.");
+                    return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]];
+                }
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::EpochSolution const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochSolution(std::map<Epoch, EpochSolution const*> const& solutions, Epoch const& epoch) {
+                    auto epochSolutionIt = solutions.find(epoch);
+                    STORM_LOG_ASSERT(epochSolutionIt != solutions.end(), "Requested unexisting solution for epoch " << epochManager.toString(epoch) << ".");
+                    return *epochSolutionIt->second;
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState) {
+                    STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution at an unexisting product state.");
+                    STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at a state for which no solution was stored.");
+                    return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]];
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch) {
+                    STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "The model has multiple initial states.");
+                    STORM_LOG_ASSERT(!epochManager.hasBottomDimension(epoch), "Tried to get the initial state result in an epoch that still contains at least one bottom dimension.");
+                    return getStateSolution(epoch, productModel->getInitialProductState(*model.getInitialStates().begin(), model.getInitialStates()));
+                }
+    
+                template<typename ValueType, bool SingleObjectiveMode>
+                typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) {
+                    STORM_LOG_ASSERT(model.getInitialStates().get(initialStateIndex), "The given model state is not an initial state.");
+                    STORM_LOG_ASSERT(!epochManager.hasBottomDimension(epoch), "Tried to get the initial state result in an epoch that still contains at least one bottom dimension.");
+                    return getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates()));
+                }
+    
+                template<typename ValueType, bool SingleObjectiveMode>
+                EpochManager const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochManager() const {
+                    return epochManager;
+                }
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                Dimension<ValueType> const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getDimension(uint64_t dim) const {
+                    return dimensions.at(dim);
+                }
+                
+                template class MultiDimensionalRewardUnfolding<double, true>;
+                template class MultiDimensionalRewardUnfolding<double, false>;
+                template class MultiDimensionalRewardUnfolding<storm::RationalNumber, true>;
+                template class MultiDimensionalRewardUnfolding<storm::RationalNumber, false>;
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
@ -0,0 +1,128 @@
+#pragma once
+
+#include <boost/optional.hpp>
+
+#include "storm/storage/BitVector.h"
+#include "storm/storage/SparseMatrix.h"
+#include "storm/modelchecker/multiobjective/Objective.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/Dimension.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/utility/vector.h"
+#include "storm/storage/memorystructure/MemoryStructure.h"
+#include "storm/utility/Stopwatch.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                template<typename ValueType, bool SingleObjectiveMode>
+                class MultiDimensionalRewardUnfolding {
+                public:
+                    
+                    typedef typename EpochManager::Epoch Epoch; // The number of reward steps that are "left" for each dimension
+                    typedef typename EpochManager::EpochClass EpochClass;
+                    
+                    typedef typename std::conditional<SingleObjectiveMode, ValueType, std::vector<ValueType>>::type SolutionType;
+    
+                    struct EpochModel {
+                        bool epochMatrixChanged;
+                        storm::storage::SparseMatrix<ValueType> epochMatrix;
+                        storm::storage::BitVector stepChoices;
+                        std::vector<SolutionType> stepSolutions;
+                        std::vector<std::vector<ValueType>> objectiveRewards;
+                        std::vector<storm::storage::BitVector> objectiveRewardFilter;
+                        storm::storage::BitVector epochInStates;
+                    };
+                    
+                    /*
+                     *
+                     * @param model The (preprocessed) model
+                     * @param objectives The (preprocessed) objectives
+                     *
+                     */
+                    MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
+                    MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp<ValueType> const& model, std::shared_ptr<storm::logic::OperatorFormula const> objectiveFormula);
+                    
+                    ~MultiDimensionalRewardUnfolding() = default;
+                    
+                    Epoch getStartEpoch();
+                    std::vector<Epoch> getEpochComputationOrder(Epoch const& startEpoch);
+                    
+                    EpochModel& setCurrentEpoch(Epoch const& epoch);
+                    
+                    /*!
+                     * Returns the precision required for the analyzis of each epoch model in order to achieve the given overall precision
+                     */
+                    ValueType getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision);
+                    
+                    void setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions);
+                    SolutionType const& getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
+                    SolutionType const& getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);
+                    
+                    EpochManager const& getEpochManager() const;
+                    Dimension<ValueType> const& getDimension(uint64_t dim) const;
+                    
+                private:
+                
+                    void setCurrentEpochClass(Epoch const& epoch);
+                    void initialize();
+                    
+                    void initializeObjectives(std::vector<Epoch>& epochSteps);
+                    void computeMaxDimensionValues();
+                    
+                    void initializeMemoryProduct(std::vector<Epoch> const& epochSteps);
+                    
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
+                    SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const;
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
+                    SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const;
+                    
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
+                    void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const;
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
+                    void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const;
+                    
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
+                    std::string solutionToString(SolutionType const& solution) const;
+                    template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
+                    std::string solutionToString(SolutionType const& solution) const;
+                    
+                    SolutionType const& getStateSolution(Epoch const& epoch, uint64_t const& productState);
+                    struct EpochSolution {
+                        uint64_t count;
+                        std::shared_ptr<std::vector<uint64_t> const> productStateToSolutionVectorMap;
+                        std::vector<SolutionType> solutions;
+                    };
+                    std::map<Epoch, EpochSolution> epochSolutions;
+                    EpochSolution const& getEpochSolution(std::map<Epoch, EpochSolution const*> const& solutions, Epoch const& epoch);
+                    SolutionType const& getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState);
+                    
+                    storm::models::sparse::Mdp<ValueType> const& model;
+                    std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> objectives;
+                    
+                    std::unique_ptr<ProductModel<ValueType>> productModel;
+                    
+                    std::vector<uint64_t> epochModelToProductChoiceMap;
+                    std::shared_ptr<std::vector<uint64_t> const> productStateToEpochModelInStateMap;
+                    std::set<Epoch> possibleEpochSteps;
+                    
+                    EpochModel epochModel;
+                    boost::optional<Epoch> currentEpoch;
+                    
+                    EpochManager epochManager;
+                    
+                    std::vector<Dimension<ValueType>> dimensions;
+                    std::vector<storm::storage::BitVector> objectiveDimensions;
+                    
+                    
+                    storm::utility::Stopwatch swInit, swFindSol, swInsertSol, swSetEpoch, swSetEpochClass, swAux1, swAux2, swAux3, swAux4;
+                    std::vector<uint64_t> epochModelSizes;
+                    uint64_t maxSolutionsStored;
+                };
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
@ -0,0 +1,605 @@
+#include "storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h"
+
+
+#include "storm/utility/macros.h"
+#include "storm/logic/Formulas.h"
+#include "storm/logic/CloneVisitor.h"
+#include "storm/storage/memorystructure/SparseModelMemoryProduct.h"
+#include "storm/storage/memorystructure/MemoryStructureBuilder.h"
+
+#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+
+#include "storm/exceptions/UnexpectedException.h"
+#include "storm/exceptions/NotSupportedException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                template<typename ValueType>
+                ProductModel<ValueType>::ProductModel(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()) {
+                    
+                    for (uint64_t dim = 0; dim < dimensions.size(); ++dim) {
+                        if (!dimensions[dim].memoryLabel) {
+                            memoryStateManager.setDimensionWithoutMemory(dim);
+                        }
+                    }
+                    
+                    storm::storage::MemoryStructure memory = computeMemoryStructure(model, objectives);
+                    assert(memoryStateManager.getMemoryStateCount() == memory.getNumberOfStates());
+                    std::vector<MemoryState> memoryStateMap = computeMemoryStateMap(memory);
+                    
+                    storm::storage::SparseModelMemoryProduct<ValueType> productBuilder(memory.product(model));
+                    
+                    setReachableProductStates(productBuilder, originalModelSteps, memoryStateMap);
+                    product = productBuilder.build()->template as<storm::models::sparse::Mdp<ValueType>>();
+                    
+                    uint64_t numModelStates = productBuilder.getOriginalModel().getNumberOfStates();
+                    MemoryState upperMemStateBound = memoryStateManager.getUpperMemoryStateBound();
+                    uint64_t numMemoryStates = memoryStateManager.getMemoryStateCount();
+                    uint64_t numProductStates = getProduct().getNumberOfStates();
+                    
+                    // Compute a mappings from product states to model/memory states and back
+                    modelMemoryToProductStateMap.resize(upperMemStateBound * numModelStates, std::numeric_limits<uint64_t>::max());
+                    productToModelStateMap.resize(numProductStates, std::numeric_limits<uint64_t>::max());
+                    productToMemoryStateMap.resize(numProductStates, std::numeric_limits<uint64_t>::max());
+                    for (uint64_t modelState = 0; modelState < numModelStates; ++modelState) {
+                        for (uint64_t memoryStateIndex = 0; memoryStateIndex < numMemoryStates; ++memoryStateIndex) {
+                            if (productBuilder.isStateReachable(modelState, memoryStateIndex)) {
+                                uint64_t productState = productBuilder.getResultState(modelState, memoryStateIndex);
+                                modelMemoryToProductStateMap[modelState * upperMemStateBound + memoryStateMap[memoryStateIndex]] = productState;
+                                productToModelStateMap[productState] = modelState;
+                                productToMemoryStateMap[productState] = memoryStateMap[memoryStateIndex];
+                            }
+                        }
+                    }
+                    
+                    // Map choice indices of the product to the state where it origins
+                    choiceToStateMap.reserve(getProduct().getNumberOfChoices());
+                    for (uint64_t productState = 0; productState < numProductStates; ++productState) {
+                        uint64_t groupSize = getProduct().getTransitionMatrix().getRowGroupSize(productState);
+                        for (uint64_t i = 0; i < groupSize; ++i) {
+                            choiceToStateMap.push_back(productState);
+                        }
+                    }
+                    
+                    // Compute the epoch steps for the product
+                    steps.resize(getProduct().getNumberOfChoices(), 0);
+                    for (uint64_t modelState = 0; modelState < numModelStates; ++modelState) {
+                        uint64_t numChoices = productBuilder.getOriginalModel().getTransitionMatrix().getRowGroupSize(modelState);
+                        uint64_t firstChoice = productBuilder.getOriginalModel().getTransitionMatrix().getRowGroupIndices()[modelState];
+                        for (uint64_t choiceOffset = 0; choiceOffset < numChoices; ++choiceOffset) {
+                            Epoch const& step  = originalModelSteps[firstChoice + choiceOffset];
+                            if (step != 0) {
+                                for (MemoryState const& memoryState : memoryStateMap) {
+                                    if (productStateExists(modelState, memoryState)) {
+                                        uint64_t productState = getProductState(modelState, memoryState);
+                                        uint64_t productChoice = getProduct().getTransitionMatrix().getRowGroupIndices()[productState] + choiceOffset;
+                                        assert(productChoice < getProduct().getTransitionMatrix().getRowGroupIndices()[productState + 1]);
+                                        steps[productChoice] = step;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                    
+                  //  getProduct().writeDotToStream(std::cout);
+                    
+                    computeReachableStatesInEpochClasses();
+                }
+                
+                
+                template<typename ValueType>
+                storm::storage::MemoryStructure ProductModel<ValueType>::computeMemoryStructure(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const {
+                    
+                    storm::modelchecker::SparsePropositionalModelChecker<storm::models::sparse::Mdp<ValueType>> mc(model);
+                    
+                    // Create a memory structure that remembers whether (sub)objectives are satisfied
+                    storm::storage::MemoryStructure memory = storm::storage::MemoryStructureBuilder<ValueType>::buildTrivialMemoryStructure(model);
+                    for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                        if (!objectives[objIndex].formula->isProbabilityOperatorFormula()) {
+                            continue;
+                        }
+                        
+                        std::vector<uint64_t> dimensionIndexMap;
+                        for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) {
+                            dimensionIndexMap.push_back(globalDimensionIndex);
+                        }
+                        
+                        bool objectiveContainsLowerBound = false;
+                        for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) {
+                            if (!dimensions[globalDimensionIndex].isUpperBounded) {
+                                objectiveContainsLowerBound = true;
+                                break;
+                            }
+                        }
+                        
+                        // collect the memory states for this objective
+                        std::vector<storm::storage::BitVector> objMemStates;
+                        storm::storage::BitVector m(dimensionIndexMap.size(), false);
+                        for (; !m.full(); m.increment()) {
+                            objMemStates.push_back(~m);
+                        }
+                        objMemStates.push_back(~m);
+                        assert(objMemStates.size() == 1ull << dimensionIndexMap.size());
+                        
+                        // build objective memory
+                        auto objMemoryBuilder = storm::storage::MemoryStructureBuilder<ValueType>(objMemStates.size(), model);
+                        
+                        // Get the set of states that for all subobjectives satisfy either the left or the right subformula
+                        storm::storage::BitVector constraintStates(model.getNumberOfStates(), true);
+                        for (auto const& dim : objectiveDimensions[objIndex]) {
+                            auto const& dimension = dimensions[dim];
+                            STORM_LOG_ASSERT(dimension.formula->isBoundedUntilFormula(), "Unexpected Formula type");
+                            constraintStates &=
+                                    (mc.check(dimension.formula->asBoundedUntilFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector() |
+                                    mc.check(dimension.formula->asBoundedUntilFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                        }
+                        
+                        // Build the transitions between the memory states
+                        for (uint64_t memState = 0; memState < objMemStates.size(); ++memState) {
+                            auto const& memStateBV = objMemStates[memState];
+                            for (uint64_t memStatePrime = 0; memStatePrime < objMemStates.size(); ++memStatePrime) {
+                                auto const& memStatePrimeBV = objMemStates[memStatePrime];
+                                if (memStatePrimeBV.isSubsetOf(memStateBV)) {
+                                    
+                                    std::shared_ptr<storm::logic::Formula const> transitionFormula = storm::logic::Formula::getTrueFormula();
+                                    for (auto const& subObjIndex : memStateBV) {
+                                        std::shared_ptr<storm::logic::Formula const> subObjFormula = dimensions[dimensionIndexMap[subObjIndex]].formula->asBoundedUntilFormula().getRightSubformula().asSharedPointer();
+                                        if (memStatePrimeBV.get(subObjIndex)) {
+                                            subObjFormula = std::make_shared<storm::logic::UnaryBooleanStateFormula>(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, subObjFormula);
+                                        }
+                                        transitionFormula = std::make_shared<storm::logic::BinaryBooleanStateFormula>(storm::logic::BinaryBooleanStateFormula::OperatorType::And, transitionFormula, subObjFormula);
+                                    }
+                                    
+                                    storm::storage::BitVector transitionStates = mc.check(*transitionFormula)->asExplicitQualitativeCheckResult().getTruthValuesVector();
+                                    if (memStatePrimeBV.empty()) {
+                                        transitionStates |= ~constraintStates;
+                                    } else {
+                                        transitionStates &= constraintStates;
+                                    }
+                                    objMemoryBuilder.setTransition(memState, memStatePrime, transitionStates);
+                                    
+                                    // Set the initial states
+                                    if (memStateBV.full()) {
+                                        storm::storage::BitVector initialTransitionStates = model.getInitialStates() & transitionStates;
+                                        // At this point we can check whether there is an initial state that already satisfies all subObjectives.
+                                        // Such a situation is not supported as we can not reduce this (easily) to an expected reward computation.
+                                        STORM_LOG_THROW(!memStatePrimeBV.empty() || initialTransitionStates.empty() || objectiveContainsLowerBound || initialTransitionStates.isDisjointFrom(constraintStates), storm::exceptions::NotSupportedException, "The objective " << *objectives[objIndex].formula << " is already satisfied in an initial state. This special case is not supported.");
+                                        for (auto const& initState : initialTransitionStates) {
+                                            objMemoryBuilder.setInitialMemoryState(initState, memStatePrime);
+                                        }
+                                    }
+                                }
+                            }
+                        }
+    
+                        // Build the memory labels
+                        for (uint64_t memState = 0; memState < objMemStates.size(); ++memState) {
+                            auto const& memStateBV = objMemStates[memState];
+                            for (auto const& subObjIndex : memStateBV) {
+                                objMemoryBuilder.setLabel(memState, dimensions[dimensionIndexMap[subObjIndex]].memoryLabel.get());
+                            }
+                        }
+                        auto objMemory = objMemoryBuilder.build();
+                        memory = memory.product(objMemory);
+                    }
+                    return memory;
+                }
+                
+                template<typename ValueType>
+                std::vector<typename ProductModel<ValueType>::MemoryState> ProductModel<ValueType>::computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const {
+                    // Compute a mapping between the different representations of memory states
+                    std::vector<MemoryState> result;
+                    result.reserve(memory.getNumberOfStates());
+                    for (uint64_t memStateIndex = 0; memStateIndex < memory.getNumberOfStates(); ++memStateIndex) {
+                        MemoryState memState = memoryStateManager.getInitialMemoryState();
+                        std::set<std::string> stateLabels = memory.getStateLabeling().getLabelsOfState(memStateIndex);
+                        for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                            if (dimensions[dim].memoryLabel) {
+                                if (stateLabels.find(dimensions[dim].memoryLabel.get()) != stateLabels.end()) {
+                                    memoryStateManager.setRelevantDimension(memState, dim, true);
+                                } else {
+                                    memoryStateManager.setRelevantDimension(memState, dim, false);
+                                }
+                            }
+                        }
+                        result.push_back(std::move(memState));
+                    }
+                    return result;
+                }
+    
+                template<typename ValueType>
+                void ProductModel<ValueType>::setReachableProductStates(storm::storage::SparseModelMemoryProduct<ValueType>& productBuilder, std::vector<Epoch> const& originalModelSteps, std::vector<MemoryState> const& memoryStateMap) const {
+    
+                    std::vector<uint64_t> inverseMemoryStateMap(memoryStateManager.getUpperMemoryStateBound(), std::numeric_limits<uint64_t>::max());
+                    for (uint64_t memStateIndex = 0; memStateIndex < memoryStateMap.size(); ++memStateIndex) {
+                        inverseMemoryStateMap[memoryStateMap[memStateIndex]] = memStateIndex;
+                    }
+                    
+                    auto const& memory = productBuilder.getMemory();
+                    auto const& model = productBuilder.getOriginalModel();
+                    auto const& modelTransitions = model.getTransitionMatrix();
+                    
+                    std::vector<storm::storage::BitVector> reachableProductStates(memoryStateManager.getUpperMemoryStateBound());
+                    for (auto const& memState : memoryStateMap) {
+                        reachableProductStates[memState] = storm::storage::BitVector(model.getNumberOfStates(), false);
+                    }
+                    
+                    // Initialize the reachable states with the initial states
+                    EpochClass initEpochClass = epochManager.getEpochClass(epochManager.getZeroEpoch());
+                    auto memStateIt = memory.getInitialMemoryStates().begin();
+                    for (auto const& initState : model.getInitialStates()) {
+                        uint64_t transformedMemoryState = transformMemoryState(memoryStateMap[*memStateIt], initEpochClass, memoryStateManager.getInitialMemoryState());
+                        reachableProductStates[transformedMemoryState].set(initState, true);
+                        ++memStateIt;
+                    }
+                    assert(memStateIt == memory.getInitialMemoryStates().end());
+    
+                    // Find the reachable epoch classes
+                    std::set<Epoch> possibleSteps(originalModelSteps.begin(), originalModelSteps.end());
+                    std::set<EpochClass, std::function<bool(EpochClass const&, EpochClass const&)>> reachableEpochClasses(std::bind(&EpochManager::epochClassOrder, &epochManager, std::placeholders::_1, std::placeholders::_2));
+                    collectReachableEpochClasses(reachableEpochClasses, possibleSteps);
+                    
+     
+                    // Iterate over all epoch classes starting from the initial one (i.e., no bottom dimension).
+                    for (auto epochClassIt = reachableEpochClasses.rbegin(); epochClassIt != reachableEpochClasses.rend(); ++epochClassIt) {
+                        auto const& epochClass = *epochClassIt;
+                        
+                        // Find the remaining set of reachable states via DFS.
+                        std::vector<std::pair<uint64_t, MemoryState>> dfsStack;
+                        for (MemoryState const& memState : memoryStateMap) {
+                            for (auto const& modelState : reachableProductStates[memState]) {
+                                dfsStack.emplace_back(modelState, memState);
+                            }
+                        }
+                    
+                        while (!dfsStack.empty()) {
+                            uint64_t currentModelState = dfsStack.back().first;
+                            MemoryState currentMemoryState = dfsStack.back().second;
+                            uint64_t currentMemoryStateIndex = inverseMemoryStateMap[currentMemoryState];
+                            dfsStack.pop_back();
+                            
+                            for (uint64_t choice = modelTransitions.getRowGroupIndices()[currentModelState]; choice != modelTransitions.getRowGroupIndices()[currentModelState + 1]; ++choice) {
+                                
+                                for (auto transitionIt = modelTransitions.getRow(choice).begin(); transitionIt < modelTransitions.getRow(choice).end(); ++transitionIt) {
+                                    
+                                    MemoryState successorMemoryState = memoryStateMap[memory.getSuccessorMemoryState(currentMemoryStateIndex, transitionIt - modelTransitions.begin())];
+                                    successorMemoryState = transformMemoryState(successorMemoryState, epochClass, currentMemoryState);
+                                    if (!reachableProductStates[successorMemoryState].get(transitionIt->getColumn())) {
+                                        reachableProductStates[successorMemoryState].set(transitionIt->getColumn(), true);
+                                        dfsStack.emplace_back(transitionIt->getColumn(), successorMemoryState);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                    
+                    for (uint64_t memStateIndex = 0; memStateIndex < memoryStateManager.getMemoryStateCount(); ++memStateIndex) {
+                        for (auto const& modelState : reachableProductStates[memoryStateMap[memStateIndex]]) {
+                            productBuilder.addReachableState(modelState, memStateIndex);
+                        }
+                    }
+                }
+    
+                template<typename ValueType>
+                storm::models::sparse::Mdp<ValueType> const& ProductModel<ValueType>::getProduct() const {
+                    return *product;
+                }
+                
+                template<typename ValueType>
+                std::vector<typename ProductModel<ValueType>::Epoch> const& ProductModel<ValueType>::getSteps() const {
+                    return steps;
+                }
+                
+                template<typename ValueType>
+                bool ProductModel<ValueType>::productStateExists(uint64_t const& modelState, MemoryState const& memoryState) const {
+                    return modelMemoryToProductStateMap[modelState * memoryStateManager.getUpperMemoryStateBound() + memoryState] < getProduct().getNumberOfStates();
+                }
+                
+                template<typename ValueType>
+                uint64_t ProductModel<ValueType>::getProductState(uint64_t const& modelState, MemoryState const& memoryState) const {
+                    STORM_LOG_ASSERT(productStateExists(modelState, memoryState), "Tried to obtain a state in the model-memory-product that does not exist");
+                    return modelMemoryToProductStateMap[modelState * memoryStateManager.getUpperMemoryStateBound() + memoryState];
+                }
+    
+                template<typename ValueType>
+                uint64_t ProductModel<ValueType>::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const {
+                    auto productInitStateIt = getProduct().getInitialStates().begin();
+                    productInitStateIt += initialModelStates.getNumberOfSetBitsBeforeIndex(initialModelState);
+                    STORM_LOG_ASSERT(getModelState(*productInitStateIt) == initialModelState, "Could not find the corresponding initial state in the product model.");
+                    return transformProductState(*productInitStateIt, epochManager.getEpochClass(epochManager.getZeroEpoch()), memoryStateManager.getInitialMemoryState());
+                }
+                
+                template<typename ValueType>
+                uint64_t ProductModel<ValueType>::getModelState(uint64_t const& productState) const {
+                    return productToModelStateMap[productState];
+                }
+                
+                template<typename ValueType>
+                typename ProductModel<ValueType>::MemoryState ProductModel<ValueType>::getMemoryState(uint64_t const& productState) const {
+                    return productToMemoryStateMap[productState];
+                }
+                
+                template<typename ValueType>
+                MemoryStateManager const& ProductModel<ValueType>::getMemoryStateManager() const {
+                    return memoryStateManager;
+                }
+                
+                template<typename ValueType>
+                uint64_t ProductModel<ValueType>::getProductStateFromChoice(uint64_t const& productChoice) const {
+                    return choiceToStateMap[productChoice];
+                }
+                
+                template<typename ValueType>
+                std::vector<std::vector<ValueType>> ProductModel<ValueType>::computeObjectiveRewards(EpochClass const& epochClass, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const {
+                    std::vector<std::vector<ValueType>> objectiveRewards;
+                    objectiveRewards.reserve(objectives.size());
+                    
+                    for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                        auto const& formula = *objectives[objIndex].formula;
+                        if (formula.isProbabilityOperatorFormula()) {
+                            storm::modelchecker::SparsePropositionalModelChecker<storm::models::sparse::Mdp<ValueType>> mc(getProduct());
+                            std::vector<uint64_t> dimensionIndexMap;
+                            for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) {
+                                dimensionIndexMap.push_back(globalDimensionIndex);
+                            }
+                            
+                            std::shared_ptr<storm::logic::Formula const> sinkStatesFormula;
+                            for (auto const& dim : objectiveDimensions[objIndex]) {
+                                auto memLabelFormula = std::make_shared<storm::logic::AtomicLabelFormula>(dimensions[dim].memoryLabel.get());
+                                if (sinkStatesFormula) {
+                                    sinkStatesFormula = std::make_shared<storm::logic::BinaryBooleanStateFormula>(storm::logic::BinaryBooleanStateFormula::OperatorType::Or, sinkStatesFormula, memLabelFormula);
+                                } else {
+                                    sinkStatesFormula = memLabelFormula;
+                                }
+                            }
+                            sinkStatesFormula = std::make_shared<storm::logic::UnaryBooleanStateFormula>(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, sinkStatesFormula);
+                            
+                            std::vector<ValueType> objRew(getProduct().getTransitionMatrix().getRowCount(), storm::utility::zero<ValueType>());
+                            storm::storage::BitVector relevantObjectives(objectiveDimensions[objIndex].getNumberOfSetBits());
+                            
+                            while (!relevantObjectives.full()) {
+                                relevantObjectives.increment();
+                                
+                                // find out whether objective reward should be earned within this epoch class
+                                bool collectRewardInEpoch = true;
+                                for (auto const& subObjIndex : relevantObjectives) {
+                                    if (dimensions[dimensionIndexMap[subObjIndex]].isUpperBounded && epochManager.isBottomDimensionEpochClass(epochClass, dimensionIndexMap[subObjIndex])) {
+                                        collectRewardInEpoch = false;
+                                        break;
+                                    }
+                                }
+                                
+                                if (collectRewardInEpoch) {
+                                    std::shared_ptr<storm::logic::Formula const> relevantStatesFormula;
+                                    std::shared_ptr<storm::logic::Formula const> goalStatesFormula =  storm::logic::CloneVisitor().clone(*sinkStatesFormula);
+                                    for (uint64_t subObjIndex = 0; subObjIndex < dimensionIndexMap.size(); ++subObjIndex) {
+                                        std::shared_ptr<storm::logic::Formula> memLabelFormula = std::make_shared<storm::logic::AtomicLabelFormula>(dimensions[dimensionIndexMap[subObjIndex]].memoryLabel.get());
+                                        if (relevantObjectives.get(subObjIndex)) {
+                                            auto rightSubFormula = dimensions[dimensionIndexMap[subObjIndex]].formula->asBoundedUntilFormula().getRightSubformula().asSharedPointer();
+                                            goalStatesFormula = std::make_shared<storm::logic::BinaryBooleanStateFormula>(storm::logic::BinaryBooleanStateFormula::OperatorType::And, goalStatesFormula, rightSubFormula);
+                                        } else {
+                                            memLabelFormula = std::make_shared<storm::logic::UnaryBooleanStateFormula>(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, memLabelFormula);
+                                        }
+                                        if (relevantStatesFormula) {
+                                            relevantStatesFormula = std::make_shared<storm::logic::BinaryBooleanStateFormula>(storm::logic::BinaryBooleanStateFormula::OperatorType::And, relevantStatesFormula, memLabelFormula);
+                                        } else {
+                                            relevantStatesFormula = memLabelFormula;
+                                        }
+                                    }
+                                    
+                                    storm::storage::BitVector relevantStates = mc.check(*relevantStatesFormula)->asExplicitQualitativeCheckResult().getTruthValuesVector();
+                                    storm::storage::BitVector relevantChoices = getProduct().getTransitionMatrix().getRowFilter(relevantStates);
+                                    storm::storage::BitVector goalStates = mc.check(*goalStatesFormula)->asExplicitQualitativeCheckResult().getTruthValuesVector();
+                                    for (auto const& choice : relevantChoices) {
+                                        objRew[choice] += getProduct().getTransitionMatrix().getConstrainedRowSum(choice, goalStates);
+                                    }
+                                }
+                            }
+                            
+                            objectiveRewards.push_back(std::move(objRew));
+                            
+                        } else if (formula.isRewardOperatorFormula()) {
+                            auto const& rewModel = getProduct().getRewardModel(formula.asRewardOperatorFormula().getRewardModelName());
+                            STORM_LOG_THROW(!rewModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Reward model has transition rewards which is not expected.");
+                            bool rewardCollectedInEpoch = true;
+                            if (formula.getSubformula().isCumulativeRewardFormula()) {
+                                assert(objectiveDimensions[objIndex].getNumberOfSetBits() == 1);
+                                rewardCollectedInEpoch = !epochManager.isBottomDimensionEpochClass(epochClass, *objectiveDimensions[objIndex].begin());
+                            } else {
+                                STORM_LOG_THROW(formula.getSubformula().isTotalRewardFormula(), storm::exceptions::UnexpectedException, "Unexpected type of formula " << formula);
+                            }
+                            if (rewardCollectedInEpoch) {
+                                objectiveRewards.push_back(rewModel.getTotalRewardVector(getProduct().getTransitionMatrix()));
+                            } else {
+                                objectiveRewards.emplace_back(getProduct().getTransitionMatrix().getRowCount(), storm::utility::zero<ValueType>());
+                            }
+                        } else {
+                            STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected type of formula " << formula);
+                        }
+                    }
+                    
+                    return objectiveRewards;
+                }
+                
+                template<typename ValueType>
+                storm::storage::BitVector const& ProductModel<ValueType>::getInStates(EpochClass const& epochClass) const {
+                    STORM_LOG_ASSERT(inStates.find(epochClass) != inStates.end(), "Could not find InStates for the given epoch class");
+                    return inStates.find(epochClass)->second;
+                }
+    
+                template<typename ValueType>
+                void ProductModel<ValueType>::computeReachableStatesInEpochClasses() {
+                    std::set<Epoch> possibleSteps(steps.begin(), steps.end());
+                    std::set<EpochClass, std::function<bool(EpochClass const&, EpochClass const&)>> reachableEpochClasses(std::bind(&EpochManager::epochClassOrder, &epochManager, std::placeholders::_1, std::placeholders::_2));
+    
+                    collectReachableEpochClasses(reachableEpochClasses, possibleSteps);
+                    
+                    for (auto epochClassIt = reachableEpochClasses.rbegin(); epochClassIt != reachableEpochClasses.rend(); ++epochClassIt) {
+                        std::vector<EpochClass> predecessors;
+                        for (auto predecessorIt = reachableEpochClasses.rbegin(); predecessorIt != epochClassIt; ++predecessorIt) {
+                            if (epochManager.isPredecessorEpochClass(*predecessorIt, *epochClassIt)) {
+                                predecessors.push_back(*predecessorIt);
+                            }
+                        }
+                        computeReachableStates(*epochClassIt, predecessors);
+                    }
+                }
+    
+                template<typename ValueType>
+                void ProductModel<ValueType>::collectReachableEpochClasses(std::set<EpochClass, std::function<bool(EpochClass const&, EpochClass const&)>>& reachableEpochClasses, std::set<Epoch> const& possibleSteps) const {
+                    
+                    Epoch startEpoch = epochManager.getZeroEpoch();
+                    for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                        STORM_LOG_ASSERT(dimensions[dim].maxValue,  "No max-value for dimension " << dim << " was given.");
+                        epochManager.setDimensionOfEpoch(startEpoch, dim, dimensions[dim].maxValue.get());
+                    }
+                    
+                    std::set<Epoch> seenEpochs({startEpoch});
+                    std::vector<Epoch> dfsStack({startEpoch});
+                    
+                    reachableEpochClasses.insert(epochManager.getEpochClass(startEpoch));
+    
+                    // Perform a DFS to find all the reachable epochs
+                    while (!dfsStack.empty()) {
+                        Epoch currentEpoch = dfsStack.back();
+                        dfsStack.pop_back();
+                        for (auto const& step : possibleSteps) {
+                            Epoch successorEpoch = epochManager.getSuccessorEpoch(currentEpoch, step);
+                            if (seenEpochs.insert(successorEpoch).second) {
+                                reachableEpochClasses.insert(epochManager.getEpochClass(successorEpoch));
+                                dfsStack.push_back(std::move(successorEpoch));
+                            }
+                        }
+                    }
+                }
+    
+                template<typename ValueType>
+                void ProductModel<ValueType>::computeReachableStates(EpochClass const& epochClass, std::vector<EpochClass> const& predecessors) {
+                    
+                    storm::storage::BitVector bottomDimensions(epochManager.getDimensionCount(), false);
+                    for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                        if (epochManager.isBottomDimensionEpochClass(epochClass, dim)) {
+                            bottomDimensions.set(dim, true);
+                        }
+                    }
+                    storm::storage::BitVector nonBottomDimensions = ~bottomDimensions;
+                    
+                    storm::storage::BitVector ecInStates(getProduct().getNumberOfStates(), false);
+                    
+                    if (!epochManager.hasBottomDimensionEpochClass(epochClass)) {
+                        for (auto const& initState : getProduct().getInitialStates()) {
+                            uint64_t transformedInitState = transformProductState(initState, epochClass, memoryStateManager.getInitialMemoryState());
+                            ecInStates.set(transformedInitState, true);
+                        }
+                    }
+                    
+                    for (auto const& predecessor : predecessors) {
+                        storm::storage::BitVector positiveStepDimensions(epochManager.getDimensionCount(), false);
+                        for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
+                            if (!epochManager.isBottomDimensionEpochClass(predecessor, dim) && bottomDimensions.get(dim)) {
+                                positiveStepDimensions.set(dim, true);
+                            }
+                        }
+                        STORM_LOG_ASSERT(reachableStates.find(predecessor) != reachableStates.end(), "Could not find reachable states of predecessor epoch class.");
+                        storm::storage::BitVector predecessorStates = reachableStates.find(predecessor)->second;
+                        for (auto const& predecessorState : predecessorStates) {
+                            uint64_t predecessorMemoryState = getMemoryState(predecessorState);
+                            for (uint64_t choice = getProduct().getTransitionMatrix().getRowGroupIndices()[predecessorState]; choice < getProduct().getTransitionMatrix().getRowGroupIndices()[predecessorState + 1]; ++choice) {
+                                bool choiceLeadsToThisClass = false;
+                                Epoch const& choiceStep = getSteps()[choice];
+                                for (auto const& dim : positiveStepDimensions) {
+                                    if (epochManager.getDimensionOfEpoch(choiceStep, dim) > 0) {
+                                        choiceLeadsToThisClass = true;
+                                    }
+                                }
+                                
+                                if (choiceLeadsToThisClass) {
+                                    for (auto const& transition : getProduct().getTransitionMatrix().getRow(choice)) {
+                                        uint64_t successorState = transformProductState(transition.getColumn(), epochClass, predecessorMemoryState);
+    
+                                        ecInStates.set(successorState, true);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                    
+                    // Find all states reachable from an InState via DFS.
+                    storm::storage::BitVector ecReachableStates = ecInStates;
+                    std::vector<uint64_t> dfsStack(ecReachableStates.begin(), ecReachableStates.end());
+                    
+                    while (!dfsStack.empty()) {
+                        uint64_t currentState = dfsStack.back();
+                        uint64_t currentMemoryState = getMemoryState(currentState);
+                        dfsStack.pop_back();
+                        
+                        for (uint64_t choice = getProduct().getTransitionMatrix().getRowGroupIndices()[currentState]; choice != getProduct().getTransitionMatrix().getRowGroupIndices()[currentState + 1]; ++choice) {
+                            
+                            bool choiceLeadsOutsideOfEpoch = false;
+                            Epoch const& choiceStep = getSteps()[choice];
+                            for (auto const& dim : nonBottomDimensions) {
+                                if (epochManager.getDimensionOfEpoch(choiceStep, dim) > 0) {
+                                    choiceLeadsOutsideOfEpoch = true;
+                                    break;
+                                }
+                            }
+                            
+                            for (auto const& transition : getProduct().getTransitionMatrix().getRow(choice)) {
+                                uint64_t successorState = transformProductState(transition.getColumn(), epochClass, currentMemoryState);
+                                if (choiceLeadsOutsideOfEpoch) {
+                                    ecInStates.set(successorState, true);
+                                }
+                                if (!ecReachableStates.get(successorState)) {
+                                    ecReachableStates.set(successorState, true);
+                                    dfsStack.push_back(successorState);
+                                }
+                            }
+                        }
+                    }
+                    
+                    reachableStates[epochClass] = std::move(ecReachableStates);
+                    
+                    inStates[epochClass] = std::move(ecInStates);
+                }
+    
+                template<typename ValueType>
+                typename ProductModel<ValueType>::MemoryState ProductModel<ValueType>::transformMemoryState(MemoryState const& memoryState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const {
+                    MemoryState memoryStatePrime = memoryState;
+                    
+                    for (auto const& objDimensions : objectiveDimensions) {
+                        for (auto const& dim : objDimensions) {
+                            auto const& dimension = dimensions[dim];
+                            if (dimension.memoryLabel) {
+                                bool dimUpperBounded = dimension.isUpperBounded;
+                                bool dimBottom = epochManager.isBottomDimensionEpochClass(epochClass, dim);
+                                if (dimUpperBounded && dimBottom && memoryStateManager.isRelevantDimension(predecessorMemoryState, dim)) {
+                                    STORM_LOG_ASSERT(objDimensions == dimension.dependentDimensions, "Unexpected set of dependent dimensions");
+                                    memoryStateManager.setRelevantDimensions(memoryStatePrime, objDimensions, false);
+                                    break;
+                                } else if (!dimUpperBounded && !dimBottom && memoryStateManager.isRelevantDimension(predecessorMemoryState, dim)) {
+                                    memoryStateManager.setRelevantDimensions(memoryStatePrime, dimension.dependentDimensions, true);
+                                }
+                            }
+                        }
+                    }
+                    
+                  //  std::cout << "Transformed memory state " << memoryStateManager.toString(memoryState) << " at epoch class " << epochClass << " with predecessor " << memoryStateManager.toString(predecessorMemoryState) << " to " << memoryStateManager.toString(memoryStatePrime) << std::endl;
+        
+                    return memoryStatePrime;
+                }
+    
+                template<typename ValueType>
+                uint64_t ProductModel<ValueType>::transformProductState(uint64_t const& productState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const {
+                    return getProductState(getModelState(productState), transformMemoryState(getMemoryState(productState), epochClass, predecessorMemoryState));
+                }
+                
+                template class ProductModel<double>;
+                template class ProductModel<storm::RationalNumber>;
+            }
+        }
+    }
+}
--- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
+++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
@ -0,0 +1,83 @@
+#pragma once
+
+#include <boost/optional.hpp>
+
+#include "storm/storage/BitVector.h"
+#include "storm/modelchecker/multiobjective/Objective.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h"
+#include "storm/modelchecker/prctl/helper/rewardbounded/Dimension.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/utility/vector.h"
+#include "storm/storage/memorystructure/MemoryStructure.h"
+#include "storm/storage/memorystructure/SparseModelMemoryProduct.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            namespace rewardbounded {
+                
+                template<typename ValueType>
+                class ProductModel {
+                public:
+                    
+                    typedef typename EpochManager::Epoch Epoch;
+                    typedef typename EpochManager::EpochClass EpochClass;
+                    typedef typename MemoryStateManager::MemoryState MemoryState;
+    
+                    
+                    ProductModel(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps);
+                    
+                    storm::models::sparse::Mdp<ValueType> const& getProduct() const;
+                    std::vector<Epoch> const& getSteps() const;
+                    
+                    bool productStateExists(uint64_t const& modelState, uint64_t const& memoryState) const;
+                    uint64_t getProductState(uint64_t const& modelState, uint64_t const& memoryState) const;
+                    uint64_t getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const;
+                    uint64_t getModelState(uint64_t const& productState) const;
+                    MemoryState getMemoryState(uint64_t const& productState) const;
+                    MemoryStateManager const& getMemoryStateManager() const;
+                    
+                    uint64_t getProductStateFromChoice(uint64_t const& productChoice) const;
+                    
+                    std::vector<std::vector<ValueType>> computeObjectiveRewards(EpochClass const& epochClass, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const;
+                    storm::storage::BitVector const& getInStates(EpochClass const& epochClass) const;
+    
+                    
+                    MemoryState transformMemoryState(MemoryState const& memoryState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
+                    uint64_t transformProductState(uint64_t const& productState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
+    
+                    
+                private:
+                    
+                    storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Mdp<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) const;
+                    std::vector<MemoryState> computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const;
+    
+                    
+                    void setReachableProductStates(storm::storage::SparseModelMemoryProduct<ValueType>& productBuilder, std::vector<Epoch> const& originalModelSteps, std::vector<MemoryState> const& memoryStateMap) const;
+                    
+                    void collectReachableEpochClasses(std::set<EpochClass, std::function<bool(EpochClass const&, EpochClass const&)>>& reachableEpochClasses, std::set<Epoch> const& possibleSteps) const;
+                    
+                    void computeReachableStatesInEpochClasses();
+                    void computeReachableStates(EpochClass const& epochClass, std::vector<EpochClass> const& predecessors);
+       
+                    std::vector<Dimension<ValueType>> const& dimensions;
+                    std::vector<storm::storage::BitVector> const& objectiveDimensions;
+                    EpochManager const& epochManager;
+                    MemoryStateManager memoryStateManager;
+    
+                    std::shared_ptr<storm::models::sparse::Mdp<ValueType>> product;
+                    std::vector<Epoch> steps;
+                    std::map<EpochClass, storm::storage::BitVector> reachableStates;
+                    std::map<EpochClass, storm::storage::BitVector> inStates;
+                    
+                    std::vector<uint64_t> modelMemoryToProductStateMap;
+                    std::vector<uint64_t> productToModelStateMap;
+                    std::vector<MemoryState> productToMemoryStateMap;
+                    std::vector<uint64_t> choiceToStateMap;
+                    
+                };
+            }
+        }
+    }
+}