From 86253fe88abc3931bca3f8134aa1345089c13749 Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 27 Oct 2017 13:24:59 +0200 Subject: [PATCH] moved multidimensional unfolding implementation from multiobjective into helper namespace --- .../pcaa/SparseMdpPcaaWeightVectorChecker.h | 2 +- ...dpRewardBoundedPcaaWeightVectorChecker.cpp | 7 +- ...eMdpRewardBoundedPcaaWeightVectorChecker.h | 10 +- .../prctl/SparseMdpPrctlModelChecker.cpp | 4 +- .../prctl/helper/SparseMdpPrctlHelper.cpp | 6 +- .../prctl/helper/SparseMdpPrctlHelper.h | 4 +- .../prctl/helper/rewardbounded/Dimension.h | 25 + .../helper/rewardbounded/EpochManager.cpp | 303 +++++++++ .../prctl/helper/rewardbounded/EpochManager.h | 63 ++ .../rewardbounded/MemoryStateManager.cpp | 90 +++ .../helper/rewardbounded/MemoryStateManager.h | 46 ++ .../MultiDimensionalRewardUnfolding.cpp | 640 ++++++++++++++++++ .../MultiDimensionalRewardUnfolding.h | 128 ++++ .../helper/rewardbounded/ProductModel.cpp | 605 +++++++++++++++++ .../prctl/helper/rewardbounded/ProductModel.h | 83 +++ 15 files changed, 1999 insertions(+), 17 deletions(-) create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp create mode 100644 src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h index 468175f71..28687776f 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h +++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h @@ -4,7 +4,7 @@ #include #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 { diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp index b0a703660..17d963813 100644 --- 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("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(solutionIt, solution.end()); @@ -93,14 +92,14 @@ namespace storm { } template - void SparseMdpRewardBoundedPcaaWeightVectorChecker::computeEpochSolution(typename MultiDimensionalRewardUnfolding::Epoch const& epoch, std::vector const& weightVector, EpochCheckingData& cachedData, ValueType const& precision) { + void SparseMdpRewardBoundedPcaaWeightVectorChecker::computeEpochSolution(typename helper::rewardbounded::MultiDimensionalRewardUnfolding::Epoch const& epoch, std::vector const& weightVector, EpochCheckingData& cachedData, ValueType const& precision) { ++numCheckedEpochs; swEpochModelBuild.start(); auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); swEpochModelBuild.stop(); swEpochModelAnalysis.start(); - std::vector::SolutionType> result; + std::vector::SolutionType> result; result.reserve(epochModel.epochInStates.getNumberOfSetBits()); uint64_t solutionSize = this->objectives.size() + 1; @@ -287,7 +286,7 @@ namespace storm { } template - void SparseMdpRewardBoundedPcaaWeightVectorChecker::updateCachedData(typename MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector, ValueType const& precision) { + void SparseMdpRewardBoundedPcaaWeightVectorChecker::updateCachedData(typename helper::rewardbounded::MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector, ValueType const& precision) { if (epochModel.epochMatrixChanged) { // Update the cached MinMaxSolver data diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h index a25cd0269..16ce8130d 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h +++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.h @@ -3,7 +3,7 @@ #include #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> xLinEq; std::unique_ptr> linEqSolver; - std::vector::SolutionType> solutions; + std::vector::SolutionType> solutions; }; - void computeEpochSolution(typename MultiDimensionalRewardUnfolding::Epoch const& epoch, std::vector const& weightVector, EpochCheckingData& cachedData, ValueType const& precision); + void computeEpochSolution(typename helper::rewardbounded::MultiDimensionalRewardUnfolding::Epoch const& epoch, std::vector const& weightVector, EpochCheckingData& cachedData, ValueType const& precision); - void updateCachedData(typename MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector, ValueType const& precision); + void updateCachedData(typename helper::rewardbounded::MultiDimensionalRewardUnfolding::EpochModel const& epochModel, EpochCheckingData& cachedData, std::vector const& weightVector, ValueType const& precision); storm::utility::Stopwatch swAll, swEpochModelBuild, swEpochModelAnalysis; uint64_t numSchedChanges, numCheckedEpochs, numChecks; - MultiDimensionalRewardUnfolding rewardUnfolding; + helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding; boost::optional> underApproxResult; boost::optional> overApproxResult; diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index c362b1cea..8712debb7 100644 --- 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(checkTask.getFormula().asSharedPointer(), opInfo); - storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); + helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(helper::SolutionType::UntilProbabilities, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } else { @@ -145,7 +145,7 @@ namespace storm { opInfo.bound = checkTask.getBound(); } auto formula = std::make_shared(checkTask.getFormula().asSharedPointer(), checkTask.getRewardModel(), opInfo); - storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); + helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(helper::SolutionType::ExpectedRewards, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } else { diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 0b687b97f..e63fb7374 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -84,7 +84,7 @@ namespace storm { template - std::vector analyzeTrivialEpochModel(OptimizationDirection dir, typename storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { + std::vector analyzeTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { // Assert that the epoch model is indeed trivial assert(epochModel.epochMatrix.getEntryCount() == 0); @@ -134,7 +134,7 @@ namespace storm { } template - std::vector analyzeNonTrivialEpochModel(OptimizationDirection dir, typename storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ValueType const& precision, SolutionType const& type) { + std::vector analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory 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 - std::map SparseMdpPrctlHelper::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::map SparseMdpPrctlHelper::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { storm::utility::Stopwatch swAll(true), swBuild, swCheck; auto initEpoch = rewardUnfolding.getStartEpoch(); auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch); diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h index 29438ea93..1990a4873 100644 --- 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 computeStepBoundedUntilProbabilities(storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint()); - static std::map computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, storm::modelchecker::multiobjective::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::map computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); static std::vector computeNextProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h new file mode 100644 index 000000000..80059ce1e --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/Dimension.h @@ -0,0 +1,25 @@ +#pragma once + +#include + +#include "storm/storage/BitVector.h" + +namespace storm { + namespace modelchecker { + namespace helper { + namespace rewardbounded { + + template + struct Dimension { + std::shared_ptr formula; + uint64_t objectiveIndex; + boost::optional memoryLabel; + bool isUpperBounded; + ValueType scalingFactor; + storm::storage::BitVector dependentDimensions; + boost::optional maxValue; + }; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp new file mode 100644 index 000000000..8de428483 --- /dev/null +++ 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 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 resultAsSet; + gatherPredecessorEpochs(resultAsSet, epoch, step); + return std::vector(resultAsSet.begin(), resultAsSet.end()); + } + + void EpochManager::gatherPredecessorEpochs(std::set& 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; + } + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h new file mode 100644 index 000000000..e4ddffb24 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h @@ -0,0 +1,63 @@ +#pragma once + +#include +#include +#include + +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 getPredecessorEpochs(Epoch const& epoch, Epoch const& step) const; + void gatherPredecessorEpochs(std::set& 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; + }; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.cpp new file mode 100644 index 000000000..1bbd69892 --- /dev/null +++ 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; + } + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h new file mode 100644 index 000000000..de75e422c --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include + +#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; + + }; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp new file mode 100644 index 000000000..ee09ba970 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -0,0 +1,640 @@ +#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" + +#include +#include +#include + +#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 + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::vector> const& objectives) : model(model), objectives(objectives) { + initialize(); + } + + template + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::shared_ptr 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 objective; + objective.formula = objectiveFormula; + objective.originalFormula = objective.formula; + objective.considersComplementaryEvent = false; + objectives.push_back(std::move(objective)); + + initialize(); + } + + + template + void MultiDimensionalRewardUnfolding::initialize() { + + maxSolutionsStored = 0; + + STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives."); + std::vector epochSteps; + initializeObjectives(epochSteps); + initializeMemoryProduct(epochSteps); + } + + template + void MultiDimensionalRewardUnfolding::initializeObjectives(std::vector& epochSteps) { + std::vector> 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 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(model.getNumberOfChoices(), 1)); + dimension.scalingFactor = storm::utility::one(); + } 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 actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix()); + auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(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 dimension; + dimension.formula = formula.getSubformula().asSharedPointer(); + dimension.objectiveIndex = objIndex; + dimension.isUpperBounded = true; + if (subformula.getTimeBoundReference().isTimeBound() || subformula.getTimeBoundReference().isStepBound()) { + dimensionWiseEpochSteps.push_back(std::vector(model.getNumberOfChoices(), 1)); + dimension.scalingFactor = storm::utility::one(); + } 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 actionRewards = rewardModel.getTotalRewardVector(this->model.getTransitionMatrix()); + auto discretizedRewardsAndFactor = storm::utility::vector::toIntegralVector(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 + void MultiDimensionalRewardUnfolding::initializeMemoryProduct(std::vector const& epochSteps) { + productModel = std::make_unique>(model, objectives, dimensions, objectiveDimensions, epochManager, epochSteps); + } + + template + void MultiDimensionalRewardUnfolding::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(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(); + } + } else { + discretizedBound = storm::utility::floor(discretizedBound); + } + uint64_t dimensionValue = storm::utility::convertNumber(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 MultiDimensionalRewardUnfolding::Epoch MultiDimensionalRewardUnfolding::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 + std::vector::Epoch> MultiDimensionalRewardUnfolding::getEpochComputationOrder(Epoch const& startEpoch) { + // Perform a DFS to find all the reachable epochs + std::vector dfsStack; + std::set> 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(collectedEpochs.begin(), collectedEpochs.end()); + } + + template + typename MultiDimensionalRewardUnfolding::EpochModel& MultiDimensionalRewardUnfolding::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 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 + void MultiDimensionalRewardUnfolding::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::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 const& productObjRew = productObjectiveRewards[objIndex]; + std::vector 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()); + } + 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 toEpochModelInStatesMap(productModel->getProduct().getNumberOfStates(), std::numeric_limits::max()); + std::vector epochModelStateToInStateMap = epochModel.epochInStates.getNumberOfSetBitsBeforeIndices(); + for (auto const& productState : productInStates) { + toEpochModelInStatesMap[productState] = epochModelStateToInStateMap[ecElimResult.oldToNewStateMapping[productState]]; + } + productStateToEpochModelInStateMap = std::make_shared 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 + template::type> + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const { + return solution * scalingFactor; + } + + template + template::type> + typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::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 + template::type> + void MultiDimensionalRewardUnfolding::addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const { + solution += solutionToAdd * scalingFactor; + } + + template + template::type> + void MultiDimensionalRewardUnfolding::addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const { + storm::utility::vector::addScaledVector(solution, solutionToAdd, scalingFactor); + } + + template + template::type> + std::string MultiDimensionalRewardUnfolding::solutionToString(SolutionType const& solution) const { + std::stringstream stringstream; + stringstream << solution; + return stringstream.str(); + } + + template + template::type> + std::string MultiDimensionalRewardUnfolding::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 + ValueType MultiDimensionalRewardUnfolding::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) { + // if (storm::settings::getModule().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(sumOfDimensions); + } + + template + void MultiDimensionalRewardUnfolding::setSolutionForCurrentEpoch(std::vector&& 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 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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::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 MultiDimensionalRewardUnfolding::EpochSolution const& MultiDimensionalRewardUnfolding::getEpochSolution(std::map 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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::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 MultiDimensionalRewardUnfolding::SolutionType const& MultiDimensionalRewardUnfolding::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 + EpochManager const& MultiDimensionalRewardUnfolding::getEpochManager() const { + return epochManager; + } + + template + Dimension const& MultiDimensionalRewardUnfolding::getDimension(uint64_t dim) const { + return dimensions.at(dim); + } + + template class MultiDimensionalRewardUnfolding; + template class MultiDimensionalRewardUnfolding; + template class MultiDimensionalRewardUnfolding; + template class MultiDimensionalRewardUnfolding; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h new file mode 100644 index 000000000..ea6c23080 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -0,0 +1,128 @@ +#pragma once + +#include + +#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 + 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>::type SolutionType; + + struct EpochModel { + bool epochMatrixChanged; + storm::storage::SparseMatrix epochMatrix; + storm::storage::BitVector stepChoices; + std::vector stepSolutions; + std::vector> objectiveRewards; + std::vector objectiveRewardFilter; + storm::storage::BitVector epochInStates; + }; + + /* + * + * @param model The (preprocessed) model + * @param objectives The (preprocessed) objectives + * + */ + MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::vector> const& objectives); + MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::shared_ptr objectiveFormula); + + ~MultiDimensionalRewardUnfolding() = default; + + Epoch getStartEpoch(); + std::vector 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&& 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 const& getDimension(uint64_t dim) const; + + private: + + void setCurrentEpochClass(Epoch const& epoch); + void initialize(); + + void initializeObjectives(std::vector& epochSteps); + void computeMaxDimensionValues(); + + void initializeMemoryProduct(std::vector const& epochSteps); + + template::type = 0> + SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const; + template::type = 0> + SolutionType getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const; + + template::type = 0> + void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const; + template::type = 0> + void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const; + + template::type = 0> + std::string solutionToString(SolutionType const& solution) const; + template::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 const> productStateToSolutionVectorMap; + std::vector solutions; + }; + std::map epochSolutions; + EpochSolution const& getEpochSolution(std::map const& solutions, Epoch const& epoch); + SolutionType const& getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState); + + storm::models::sparse::Mdp const& model; + std::vector> objectives; + + std::unique_ptr> productModel; + + std::vector epochModelToProductChoiceMap; + std::shared_ptr const> productStateToEpochModelInStateMap; + std::set possibleEpochSteps; + + EpochModel epochModel; + boost::optional currentEpoch; + + EpochManager epochManager; + + std::vector> dimensions; + std::vector objectiveDimensions; + + + storm::utility::Stopwatch swInit, swFindSol, swInsertSol, swSetEpoch, swSetEpochClass, swAux1, swAux2, swAux3, swAux4; + std::vector epochModelSizes; + uint64_t maxSolutionsStored; + }; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp new file mode 100644 index 000000000..8b4223fd2 --- /dev/null +++ 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 + ProductModel::ProductModel(storm::models::sparse::Mdp const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector 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 memoryStateMap = computeMemoryStateMap(memory); + + storm::storage::SparseModelMemoryProduct productBuilder(memory.product(model)); + + setReachableProductStates(productBuilder, originalModelSteps, memoryStateMap); + product = productBuilder.build()->template as>(); + + 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::max()); + productToModelStateMap.resize(numProductStates, std::numeric_limits::max()); + productToMemoryStateMap.resize(numProductStates, std::numeric_limits::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 + storm::storage::MemoryStructure ProductModel::computeMemoryStructure(storm::models::sparse::Mdp const& model, std::vector> const& objectives) const { + + storm::modelchecker::SparsePropositionalModelChecker> mc(model); + + // Create a memory structure that remembers whether (sub)objectives are satisfied + storm::storage::MemoryStructure memory = storm::storage::MemoryStructureBuilder::buildTrivialMemoryStructure(model); + for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { + if (!objectives[objIndex].formula->isProbabilityOperatorFormula()) { + continue; + } + + std::vector 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 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(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 transitionFormula = storm::logic::Formula::getTrueFormula(); + for (auto const& subObjIndex : memStateBV) { + std::shared_ptr subObjFormula = dimensions[dimensionIndexMap[subObjIndex]].formula->asBoundedUntilFormula().getRightSubformula().asSharedPointer(); + if (memStatePrimeBV.get(subObjIndex)) { + subObjFormula = std::make_shared(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, subObjFormula); + } + transitionFormula = std::make_shared(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 + std::vector::MemoryState> ProductModel::computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const { + // Compute a mapping between the different representations of memory states + std::vector result; + result.reserve(memory.getNumberOfStates()); + for (uint64_t memStateIndex = 0; memStateIndex < memory.getNumberOfStates(); ++memStateIndex) { + MemoryState memState = memoryStateManager.getInitialMemoryState(); + std::set 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 + void ProductModel::setReachableProductStates(storm::storage::SparseModelMemoryProduct& productBuilder, std::vector const& originalModelSteps, std::vector const& memoryStateMap) const { + + std::vector inverseMemoryStateMap(memoryStateManager.getUpperMemoryStateBound(), std::numeric_limits::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 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 possibleSteps(originalModelSteps.begin(), originalModelSteps.end()); + std::set> 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> 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 + storm::models::sparse::Mdp const& ProductModel::getProduct() const { + return *product; + } + + template + std::vector::Epoch> const& ProductModel::getSteps() const { + return steps; + } + + template + bool ProductModel::productStateExists(uint64_t const& modelState, MemoryState const& memoryState) const { + return modelMemoryToProductStateMap[modelState * memoryStateManager.getUpperMemoryStateBound() + memoryState] < getProduct().getNumberOfStates(); + } + + template + uint64_t ProductModel::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 + uint64_t ProductModel::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 + uint64_t ProductModel::getModelState(uint64_t const& productState) const { + return productToModelStateMap[productState]; + } + + template + typename ProductModel::MemoryState ProductModel::getMemoryState(uint64_t const& productState) const { + return productToMemoryStateMap[productState]; + } + + template + MemoryStateManager const& ProductModel::getMemoryStateManager() const { + return memoryStateManager; + } + + template + uint64_t ProductModel::getProductStateFromChoice(uint64_t const& productChoice) const { + return choiceToStateMap[productChoice]; + } + + template + std::vector> ProductModel::computeObjectiveRewards(EpochClass const& epochClass, std::vector> const& objectives) const { + std::vector> 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> mc(getProduct()); + std::vector dimensionIndexMap; + for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) { + dimensionIndexMap.push_back(globalDimensionIndex); + } + + std::shared_ptr sinkStatesFormula; + for (auto const& dim : objectiveDimensions[objIndex]) { + auto memLabelFormula = std::make_shared(dimensions[dim].memoryLabel.get()); + if (sinkStatesFormula) { + sinkStatesFormula = std::make_shared(storm::logic::BinaryBooleanStateFormula::OperatorType::Or, sinkStatesFormula, memLabelFormula); + } else { + sinkStatesFormula = memLabelFormula; + } + } + sinkStatesFormula = std::make_shared(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, sinkStatesFormula); + + std::vector objRew(getProduct().getTransitionMatrix().getRowCount(), storm::utility::zero()); + 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 relevantStatesFormula; + std::shared_ptr goalStatesFormula = storm::logic::CloneVisitor().clone(*sinkStatesFormula); + for (uint64_t subObjIndex = 0; subObjIndex < dimensionIndexMap.size(); ++subObjIndex) { + std::shared_ptr memLabelFormula = std::make_shared(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::OperatorType::And, goalStatesFormula, rightSubFormula); + } else { + memLabelFormula = std::make_shared(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, memLabelFormula); + } + if (relevantStatesFormula) { + relevantStatesFormula = std::make_shared(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()); + } + } else { + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected type of formula " << formula); + } + } + + return objectiveRewards; + } + + template + storm::storage::BitVector const& ProductModel::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 + void ProductModel::computeReachableStatesInEpochClasses() { + std::set possibleSteps(steps.begin(), steps.end()); + std::set> 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 predecessors; + for (auto predecessorIt = reachableEpochClasses.rbegin(); predecessorIt != epochClassIt; ++predecessorIt) { + if (epochManager.isPredecessorEpochClass(*predecessorIt, *epochClassIt)) { + predecessors.push_back(*predecessorIt); + } + } + computeReachableStates(*epochClassIt, predecessors); + } + } + + template + void ProductModel::collectReachableEpochClasses(std::set>& reachableEpochClasses, std::set 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 seenEpochs({startEpoch}); + std::vector 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 + void ProductModel::computeReachableStates(EpochClass const& epochClass, std::vector 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 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 ProductModel::MemoryState ProductModel::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 + uint64_t ProductModel::transformProductState(uint64_t const& productState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const { + return getProductState(getModelState(productState), transformMemoryState(getMemoryState(productState), epochClass, predecessorMemoryState)); + } + + template class ProductModel; + template class ProductModel; + } + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h new file mode 100644 index 000000000..4aaf90796 --- /dev/null +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h @@ -0,0 +1,83 @@ +#pragma once + +#include + +#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 + class ProductModel { + public: + + typedef typename EpochManager::Epoch Epoch; + typedef typename EpochManager::EpochClass EpochClass; + typedef typename MemoryStateManager::MemoryState MemoryState; + + + ProductModel(storm::models::sparse::Mdp const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps); + + storm::models::sparse::Mdp const& getProduct() const; + std::vector 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> computeObjectiveRewards(EpochClass const& epochClass, std::vector> 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 const& model, std::vector> const& objectives) const; + std::vector computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const; + + + void setReachableProductStates(storm::storage::SparseModelMemoryProduct& productBuilder, std::vector const& originalModelSteps, std::vector const& memoryStateMap) const; + + void collectReachableEpochClasses(std::set>& reachableEpochClasses, std::set const& possibleSteps) const; + + void computeReachableStatesInEpochClasses(); + void computeReachableStates(EpochClass const& epochClass, std::vector const& predecessors); + + std::vector> const& dimensions; + std::vector const& objectiveDimensions; + EpochManager const& epochManager; + MemoryStateManager memoryStateManager; + + std::shared_ptr> product; + std::vector steps; + std::map reachableStates; + std::map inStates; + + std::vector modelMemoryToProductStateMap; + std::vector productToModelStateMap; + std::vector productToMemoryStateMap; + std::vector choiceToStateMap; + + }; + } + } + } +} \ No newline at end of file