diff --git a/src/storm-pomdp-cli/settings/modules/POMDPSettings.cpp b/src/storm-pomdp-cli/settings/modules/POMDPSettings.cpp index 617cebf36..08f1e8d50 100644 --- a/src/storm-pomdp-cli/settings/modules/POMDPSettings.cpp +++ b/src/storm-pomdp-cli/settings/modules/POMDPSettings.cpp @@ -15,6 +15,7 @@ namespace storm { const std::string POMDPSettings::moduleName = "pomdp"; const std::string exportAsParametricModelOption = "parametric-drn"; const std::string gridApproximationOption = "gridapproximation"; + const std::string limitBeliefExplorationOption = "limit-exploration"; const std::string qualitativeReductionOption = "qualitativereduction"; const std::string analyzeUniqueObservationsOption = "uniqueobservations"; const std::string mecReductionOption = "mecreduction"; @@ -38,13 +39,9 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, fscmode, false, "Sets the way the pMC is obtained").addArgument(storm::settings::ArgumentBuilder::createStringArgument("type", "type name").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(fscModes)).setDefaultValueString("standard").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, transformBinaryOption, false, "Transforms the pomdp to a binary pomdp.").build()); this->addOption(storm::settings::OptionBuilder(moduleName, transformSimpleOption, false, "Transforms the pomdp to a binary and simple pomdp.").build()); - this->addOption(storm::settings::OptionBuilder(moduleName, gridApproximationOption, false, - "Analyze the POMDP using grid approximation.").addArgument( - storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("resolution", - "the resolution of the grid").setDefaultValueUnsignedInteger( - 10).addValidatorUnsignedInteger( - storm::settings::ArgumentValidatorFactory::createUnsignedGreaterValidator( - 0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, gridApproximationOption, false,"Analyze the POMDP using grid approximation.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("resolution","the resolution of the grid").setDefaultValueUnsignedInteger(10).addValidatorUnsignedInteger(storm::settings::ArgumentValidatorFactory::createUnsignedGreaterValidator(0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, limitBeliefExplorationOption, false,"Sets whether to early in the belief space exploration if upper and lower bound are close").addArgument( + storm::settings::ArgumentBuilder::createDoubleArgument("threshold","the difference between upper and lower bound when to stop").setDefaultValueDouble(0.0).addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0,1)).build()).build()); } bool POMDPSettings::isExportToParametricSet() const { @@ -76,10 +73,17 @@ namespace storm { } uint64_t POMDPSettings::getGridResolution() const { - return this->getOption(gridApproximationOption).getArgumentByName( - "resolution").getValueAsUnsignedInteger(); + return this->getOption(gridApproximationOption).getArgumentByName("resolution").getValueAsUnsignedInteger(); } - + + bool POMDPSettings::isLimitExplorationSet() const { + return this->getOption(limitBeliefExplorationOption).getHasOptionBeenSet(); + } + + double POMDPSettings::getExplorationThreshold() const { + return this->getOption(limitBeliefExplorationOption).getArgumentByName("threshold").getValueAsDouble(); + } + uint64_t POMDPSettings::getMemoryBound() const { return this->getOption(memoryBoundOption).getArgumentByName("bound").getValueAsUnsignedInteger(); } diff --git a/src/storm-pomdp-cli/settings/modules/POMDPSettings.h b/src/storm-pomdp-cli/settings/modules/POMDPSettings.h index 1e68871e8..c0c9ecf44 100644 --- a/src/storm-pomdp-cli/settings/modules/POMDPSettings.h +++ b/src/storm-pomdp-cli/settings/modules/POMDPSettings.h @@ -27,6 +27,7 @@ namespace storm { bool isQualitativeReductionSet() const; bool isGridApproximationSet() const; + bool isLimitExplorationSet() const; bool isAnalyzeUniqueObservationsSet() const; bool isMecReductionSet() const; bool isSelfloopReductionSet() const; @@ -36,6 +37,7 @@ namespace storm { uint64_t getMemoryBound() const; uint64_t getGridResolution() const; + double getExplorationThreshold() const; storm::storage::PomdpMemoryPattern getMemoryPattern() const; bool check() const override; diff --git a/src/storm-pomdp-cli/storm-pomdp.cpp b/src/storm-pomdp-cli/storm-pomdp.cpp index 0af965ddd..8b28da0bf 100644 --- a/src/storm-pomdp-cli/storm-pomdp.cpp +++ b/src/storm-pomdp-cli/storm-pomdp.cpp @@ -203,7 +203,7 @@ int main(const int argc, const char** argv) { //result = checker.refineReachabilityProbability(*pomdp, targetObservationSet,probFormula.getOptimalityType() == storm::OptimizationDirection::Minimize, pomdpSettings.getGridResolution(),1,10); result = checker.computeReachabilityProbabilityOTF(*pomdp, targetObservationSet, probFormula.getOptimalityType() == storm::OptimizationDirection::Minimize, - pomdpSettings.getGridResolution()); + pomdpSettings.getGridResolution(), pomdpSettings.getExplorationThreshold()); overRes = result->OverapproximationValue; underRes = result->UnderapproximationValue; if (overRes != underRes) { diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp index 3939913c6..fd8502563 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp @@ -52,20 +52,9 @@ namespace storm { std::unique_ptr> ApproximatePOMDPModelchecker::computeReachabilityOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, uint64_t gridResolution, - bool computeRewards) { + bool computeRewards, double explorationThreshold) { //TODO For the prototypical implementation, I put the refinement loop here. I'll change this later on - - // we define some positional scheduler for the POMDP as an experimental lower bound - storm::storage::Scheduler pomdpScheduler(pomdp.getNumberOfStates()); - for (uint32_t obs = 0; obs < pomdp.getNrObservations(); ++obs) { - auto obsStates = pomdp.getStatesWithObservation(obs); - // select a random action for all states with the same observation - uint64_t chosenAction = std::rand() % pomdp.getNumberOfChoices(obsStates.front()); - for (auto const &state : obsStates) { - pomdpScheduler.setChoice(chosenAction, state); - } - } - + storm::utility::Stopwatch underlyingWatch(true); storm::models::sparse::StateLabeling underlyingMdpLabeling(pomdp.getStateLabeling()); underlyingMdpLabeling.addLabel("goal"); std::vector goalStates; @@ -84,13 +73,23 @@ namespace storm { std::shared_ptr underlyingProperty = storm::api::extractFormulasFromProperties(propVector).front(); STORM_PRINT("Underlying MDP" << std::endl) underlyingMdp.printModelInformationToStream(std::cout); - - std::unique_ptr underlyingRes( - storm::api::verifyWithSparseEngine(underlyingModel, storm::api::createTask(underlyingProperty, false))); + std::unique_ptr underlyingRes(storm::api::verifyWithSparseEngine(underlyingModel, storm::api::createTask(underlyingProperty, false))); STORM_LOG_ASSERT(underlyingRes, "Result not exist."); underlyingRes->filter(storm::modelchecker::ExplicitQualitativeCheckResult(storm::storage::BitVector(underlyingMdp.getNumberOfStates(), true))); auto mdpResultMap = underlyingRes->asExplicitQuantitativeCheckResult().getValueMap(); + underlyingWatch.stop(); + storm::utility::Stopwatch positionalWatch(true); + // we define some positional scheduler for the POMDP as an experimental lower bound + storm::storage::Scheduler pomdpScheduler(pomdp.getNumberOfStates()); + for (uint32_t obs = 0; obs < pomdp.getNrObservations(); ++obs) { + auto obsStates = pomdp.getStatesWithObservation(obs); + // select a random action for all states with the same observation + uint64_t chosenAction = std::rand() % pomdp.getNumberOfChoices(obsStates.front()); + for (auto const &state : obsStates) { + pomdpScheduler.setChoice(chosenAction, state); + } + } auto underApproxModel = underlyingMdp.applyScheduler(pomdpScheduler, false); STORM_PRINT("Random Positional Scheduler" << std::endl) underApproxModel->printModelInformationToStream(std::cout); @@ -99,6 +98,9 @@ namespace storm { STORM_LOG_ASSERT(underapproxRes, "Result not exist."); underapproxRes->filter(storm::modelchecker::ExplicitQualitativeCheckResult(storm::storage::BitVector(underApproxModel->getNumberOfStates(), true))); auto mdpUnderapproxResultMap = underapproxRes->asExplicitQuantitativeCheckResult().getValueMap(); + positionalWatch.stop(); + + STORM_PRINT("Preprocessing Times: " << underlyingWatch << " / " << positionalWatch << std::endl) STORM_PRINT("Use On-The-Fly Grid Generation" << std::endl) std::vector> beliefList; @@ -118,25 +120,30 @@ namespace storm { std::map>> nextBelieves; // current ID -> action -> reward std::map> beliefActionRewards; - - std::vector hintVector; uint64_t nextId = 0; storm::utility::Stopwatch expansionTimer(true); - // Initial belief always has ID 0 + // Initial belief always has belief ID 0 storm::pomdp::Belief initialBelief = getInitialBelief(pomdp, nextId); ++nextId; beliefList.push_back(initialBelief); beliefIsTarget.push_back(targetObservations.find(initialBelief.observation) != targetObservations.end()); - // These are the components to build the MDP from the grid + // These are the components to build the MDP from the grid TODO make a better stucture to allow for fast reverse lookups (state-->belief) as it is a bijective function (boost:bimap?) std::map beliefStateMap; - std::vector>> mdpTransitions; - std::vector targetStates; - uint64_t mdpStateId = 0; + // Reserve states 0 and 1 as always sink/goal states + std::vector>> mdpTransitions = {{{{0, storm::utility::one()}}},{{{1, storm::utility::one()}}}}; + // Hint vector for the MDP modelchecker (initialize with constant sink/goal values) + std::vector hintVector = {storm::utility::zero(), storm::utility::one()}; + std::vector targetStates = {1}; + uint64_t mdpStateId = 2; beliefStateMap[initialBelief.id] = mdpStateId; ++mdpStateId; + // Map to save the weighted values resulting from the preprocessing for the beliefs / indices in beliefSpace + std::map weightedSumOverMap; + std::map weightedSumUnderMap; + // for the initial belief, add the triangulated initial states std::pair>, std::vector> initTemp = computeSubSimplexAndLambdas(initialBelief.probabilities, gridResolution); std::vector> initSubSimplex = initTemp.first; @@ -145,8 +152,6 @@ namespace storm { subSimplexCache[0] = initSubSimplex; lambdaCache[0] = initLambdas; } - - std::vector> initTransitionsInBelief; std::map initTransitionInActionBelief; bool initInserted = false; @@ -156,6 +161,15 @@ namespace storm { if (searchResult == uint64_t(-1) || (searchResult == 0 && !initInserted)) { if (searchResult == 0) { // the initial belief is on the grid itself + auto tempWeightedSumOver = storm::utility::zero(); + auto tempWeightedSumUnder = storm::utility::zero(); + for (uint64_t i = 0; i < initSubSimplex[j].size(); ++i) { + tempWeightedSumOver += initSubSimplex[j][i] * storm::utility::convertNumber(mdpResultMap[i]); + tempWeightedSumUnder += initSubSimplex[j][i] * storm::utility::convertNumber(mdpUnderapproxResultMap[i]); + } + weightedSumOverMap[initialBelief.id] = tempWeightedSumOver; + weightedSumUnderMap[initialBelief.id] = tempWeightedSumUnder; + initInserted = true; beliefGrid.push_back(initialBelief); beliefsToBeExpanded.push_back(0); @@ -163,6 +177,17 @@ namespace storm { : storm::utility::zero()); } else { // if the triangulated belief was not found in the list, we place it in the grid and add it to the work list + + auto tempWeightedSumOver = storm::utility::zero(); + auto tempWeightedSumUnder = storm::utility::zero(); + for (uint64_t i = 0; i < initSubSimplex[j].size(); ++i) { + tempWeightedSumOver += initSubSimplex[j][i] * storm::utility::convertNumber(mdpResultMap[i]); + tempWeightedSumUnder += initSubSimplex[j][i] * storm::utility::convertNumber(mdpUnderapproxResultMap[i]); + } + + weightedSumOverMap[nextId] = tempWeightedSumOver; + weightedSumUnderMap[nextId] = tempWeightedSumUnder; + storm::pomdp::Belief gridBelief = {nextId, initialBelief.observation, initSubSimplex[j]}; beliefList.push_back(gridBelief); beliefGrid.push_back(gridBelief); @@ -190,16 +215,22 @@ namespace storm { //beliefsToBeExpanded.push_back(initialBelief.id); I'm curious what happens if we do this instead of first triangulating. Should do nothing special if belief is on grid, otherwise it gets interesting - std::map weightedSumOverMap; - std::map weightedSumUnderMap; - - // Expand the beliefs to generate the grid on-the-fly to avoid unreachable grid points + // Expand the beliefs to generate the grid on-the-fly + if(explorationThreshold > 0){ + STORM_PRINT("Exploration threshold: " << explorationThreshold << std::endl) + } while (!beliefsToBeExpanded.empty()) { // TODO add direct generation of transition matrix uint64_t currId = beliefsToBeExpanded.front(); beliefsToBeExpanded.pop_front(); bool isTarget = beliefIsTarget[currId]; + if(cc.isLess(weightedSumOverMap[currId] - weightedSumUnderMap[currId], storm::utility::convertNumber(explorationThreshold))){ + result.emplace(std::make_pair(currId, computeRewards ? storm::utility::zero() : weightedSumOverMap[currId])); + mdpTransitions.push_back({{{1, weightedSumOverMap[currId]},{0, storm::utility::one() - weightedSumOverMap[currId]}}}); + continue; + } + if (isTarget) { // Depending on whether we compute rewards, we select the right initial result result.emplace(std::make_pair(currId, computeRewards ? storm::utility::zero() : storm::utility::one())); @@ -240,8 +271,7 @@ namespace storm { subSimplex = subSimplexCache[idNextBelief]; lambdas = lambdaCache[idNextBelief]; } else { - std::pair>, std::vector> temp = computeSubSimplexAndLambdas( - beliefList[idNextBelief].probabilities, gridResolution); + std::pair>, std::vector> temp = computeSubSimplexAndLambdas(beliefList[idNextBelief].probabilities, gridResolution); subSimplex = temp.first; lambdas = temp.second; if(cacheSubsimplices){ @@ -268,6 +298,8 @@ namespace storm { if (cc.isEqual(tempWeightedSumOver, tempWeightedSumUnder)) { hintVector.push_back(tempWeightedSumOver); + } else { + hintVector.push_back(storm::utility::zero()); } beliefsToBeExpanded.push_back(nextId); @@ -319,11 +351,10 @@ namespace storm { storm::models::sparse::StateLabeling mdpLabeling(mdpTransitions.size()); mdpLabeling.addLabel("init"); mdpLabeling.addLabel("target"); - mdpLabeling.addLabelToState("init", 0); + mdpLabeling.addLabelToState("init", beliefStateMap[initialBelief.id]); for (auto targetState : targetStates) { mdpLabeling.addLabelToState("target", targetState); } - storm::storage::sparse::ModelComponents modelComponents(buildTransitionMatrix(mdpTransitions), mdpLabeling); storm::models::sparse::Mdp overApproxMdp(modelComponents); if (computeRewards) { @@ -353,7 +384,7 @@ namespace storm { propertyString += "=? [F \"target\"]"; std::vector propertyVector = storm::api::parseProperties(propertyString); std::shared_ptr property = storm::api::extractFormulasFromProperties(propertyVector).front(); - auto task = storm::api::createTask(property, true); + auto task = storm::api::createTask(property, false); auto hint = storm::modelchecker::ExplicitModelCheckerHint(); hint.setResultHint(hintVector); auto hintPtr = std::make_shared>(hint); @@ -362,8 +393,9 @@ namespace storm { std::unique_ptr res(storm::api::verifyWithSparseEngine(model, task)); overApproxTimer.stop(); STORM_LOG_ASSERT(res, "Result not exist."); - res->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates())); - auto overApprox = res->asExplicitQuantitativeCheckResult().getValueMap().begin()->second; + res->filter(storm::modelchecker::ExplicitQualitativeCheckResult(storm::storage::BitVector(overApproxMdp.getNumberOfStates(), true))); + auto resultMap = res->asExplicitQuantitativeCheckResult().getValueMap(); + auto overApprox = resultMap[beliefStateMap[initialBelief.id]]; /* storm::utility::Stopwatch underApproxTimer(true); ValueType underApprox = computeUnderapproximationWithMDP(pomdp, beliefList, beliefIsTarget, targetObservations, observationProbabilities, nextBelieves, result, chosenActions, gridResolution, initialBelief.id, min, computeRewards); @@ -374,6 +406,11 @@ namespace storm { STORM_PRINT("Over-Approximation Result: " << overApprox << std::endl); STORM_PRINT("Under-Approximation Result: " << underApprox << std::endl); + std::map differences; + for(auto const &entry : weightedSumUnderMap){ + differences[beliefStateMap[entry.first]] = resultMap[beliefStateMap[entry.first]] - weightedSumUnderMap[entry.first]; + } + return std::make_unique>(POMDPCheckResult{overApprox, underApprox}); } @@ -499,15 +536,15 @@ namespace storm { ApproximatePOMDPModelchecker::computeReachabilityRewardOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, uint64_t gridResolution) { - return computeReachabilityOTF(pomdp, targetObservations, min, gridResolution, true); + return computeReachabilityOTF(pomdp, targetObservations, min, gridResolution, true, 0); } template std::unique_ptr> ApproximatePOMDPModelchecker::computeReachabilityProbabilityOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, - uint64_t gridResolution) { - return computeReachabilityOTF(pomdp, targetObservations, min, gridResolution, false); + uint64_t gridResolution, double explorationThreshold) { + return computeReachabilityOTF(pomdp, targetObservations, min, gridResolution, false, explorationThreshold); } template diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h index e6f935a72..438354eca 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h @@ -28,7 +28,7 @@ namespace storm { std::unique_ptr> computeReachabilityProbabilityOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, - uint64_t gridResolution); + uint64_t gridResolution, double explorationThreshold); std::unique_ptr> computeReachabilityRewardOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, @@ -57,7 +57,7 @@ namespace storm { std::unique_ptr> computeReachabilityOTF(storm::models::sparse::Pomdp const &pomdp, std::set const &targetObservations, bool min, - uint64_t gridResolution, bool computeRewards); + uint64_t gridResolution, bool computeRewards, double explorationThreshold); /** *