diff --git a/src/storm-pomdp-cli/storm-pomdp.cpp b/src/storm-pomdp-cli/storm-pomdp.cpp index 880dc29ad..8d33fc2d2 100644 --- a/src/storm-pomdp-cli/storm-pomdp.cpp +++ b/src/storm-pomdp-cli/storm-pomdp.cpp @@ -23,6 +23,7 @@ #include "storm/settings/modules/AbstractionSettings.h" #include "storm/settings/modules/BuildSettings.h" #include "storm/settings/modules/JitBuilderSettings.h" +#include "storm/settings/modules/TopologicalEquationSolverSettings.h" #include "storm/settings/modules/MultiObjectiveSettings.h" #include "storm-pomdp-cli/settings/modules/POMDPSettings.h" @@ -63,7 +64,7 @@ void initializeSettings() { storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); - + storm::settings::addModule(); storm::settings::addModule(); diff --git a/src/storm-pomdp/CMakeLists.txt b/src/storm-pomdp/CMakeLists.txt index 53643b738..37da8b794 100644 --- a/src/storm-pomdp/CMakeLists.txt +++ b/src/storm-pomdp/CMakeLists.txt @@ -17,7 +17,7 @@ set_target_properties(storm-pomdp PROPERTIES DEFINE_SYMBOL "") list(APPEND STORM_TARGETS storm-pomdp) set(STORM_TARGETS ${STORM_TARGETS} PARENT_SCOPE) -target_link_libraries(storm-pomdp PUBLIC storm ${STORM_POMDP_LINK_LIBRARIES}) +target_link_libraries(storm-pomdp PUBLIC storm storm-parsers ${STORM_POMDP_LINK_LIBRARIES}) # Install storm headers to include directory. foreach(HEADER ${STORM_POMDP_HEADERS}) diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp index e44102656..24d312c7f 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp @@ -1,6 +1,17 @@ -#include #include "ApproximatePOMDPModelchecker.h" + +#include + + +#include "storm/utility/ConstantsComparator.h" +#include "storm/models/sparse/Dtmc.h" +#include "storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "storm/utility/vector.h" +#include "storm/modelchecker/results/CheckResult.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/api/properties.h" +#include "storm-parsers/api/storm-parsers.h" namespace storm { namespace pomdp { @@ -14,35 +25,42 @@ namespace storm { /*std::unique_ptr*/ void ApproximatePOMDPModelchecker::computeReachabilityProbability( storm::models::sparse::Pomdp const &pomdp, - std::set target_observations, bool min, uint64_t gridResolution) { + std::set targetObservations, bool min, uint64_t gridResolution) { //TODO add timing uint64_t maxIterations = 100; bool finished = false; uint64_t iteration = 0; std::vector> beliefList; + std::vector beliefIsTarget; uint64_t nextId = 0; // Initial belief always has ID 0 storm::pomdp::Belief initialBelief = getInitialBelief(pomdp, nextId); ++nextId; beliefList.push_back(initialBelief); + beliefIsTarget.push_back( + targetObservations.find(initialBelief.observation) != targetObservations.end()); + std::vector> beliefGrid; - std::vector beliefIsKnown; - constructBeliefGrid(pomdp, target_observations, gridResolution, beliefList, beliefGrid, beliefIsKnown, + constructBeliefGrid(pomdp, targetObservations, gridResolution, beliefList, beliefGrid, beliefIsTarget, nextId); nextId = beliefList.size(); + // ID -> Value std::map result; std::map result_backup; + // ID -> ActionIndex std::map chosenActions; - std::vector>> observationProbabilities; - std::vector>> nextBelieves; + // ID -> Observation -> Probability + std::map>> observationProbabilities; + // current ID -> action -> next ID + std::map>> nextBelieves; for (size_t i = 0; i < beliefGrid.size(); ++i) { auto currentBelief = beliefGrid[i]; - bool isTarget = beliefIsKnown[i]; + bool isTarget = beliefIsTarget[currentBelief.id]; if (isTarget) { result.emplace(std::make_pair(currentBelief.id, storm::utility::one())); result_backup.emplace(std::make_pair(currentBelief.id, storm::utility::one())); @@ -64,6 +82,8 @@ namespace storm { uint32_t observation = iter->first; actionObservationBelieves[observation] = getBeliefAfterActionAndObservation(pomdp, beliefList, + beliefIsTarget, + targetObservations, currentBelief, action, observation, @@ -73,8 +93,9 @@ namespace storm { observationProbabilitiesInAction.push_back(actionObservationProbabilities); nextBelievesInAction.push_back(actionObservationBelieves); } - observationProbabilities.push_back(observationProbabilitiesInAction); - nextBelieves.push_back(nextBelievesInAction); + observationProbabilities.emplace( + std::make_pair(currentBelief.id, observationProbabilitiesInAction)); + nextBelieves.emplace(std::make_pair(currentBelief.id, nextBelievesInAction)); } } // Value Iteration @@ -83,9 +104,9 @@ namespace storm { STORM_LOG_DEBUG("Iteration " << iteration + 1); bool improvement = false; for (size_t i = 0; i < beliefGrid.size(); ++i) { - bool isTarget = beliefIsKnown[i]; + storm::pomdp::Belief currentBelief = beliefGrid[i]; + bool isTarget = beliefIsTarget[currentBelief.id]; if (!isTarget) { - storm::pomdp::Belief currentBelief = beliefGrid[i]; // we can take any state with the observation as they have the same number of choices uint64_t numChoices = pomdp.getNumberOfChoices( pomdp.getStatesWithObservation(currentBelief.observation).front()); @@ -97,17 +118,17 @@ namespace storm { for (uint64_t action = 0; action < numChoices; ++action) { currentValue = storm::utility::zero(); // simply change this for rewards? - for (auto iter = observationProbabilities[i][action].begin(); - iter != observationProbabilities[i][action].end(); ++iter) { + for (auto iter = observationProbabilities[currentBelief.id][action].begin(); + iter != observationProbabilities[currentBelief.id][action].end(); ++iter) { uint32_t observation = iter->first; - storm::pomdp::Belief nextBelief = beliefList[nextBelieves[i][action][observation]]; + storm::pomdp::Belief nextBelief = beliefList[nextBelieves[currentBelief.id][action][observation]]; // compute subsimplex and lambdas according to the Lovejoy paper to approximate the next belief std::pair>, std::vector> temp = computeSubSimplexAndLambdas( nextBelief.probabilities, gridResolution); std::vector> subSimplex = temp.first; std::vector lambdas = temp.second; - ValueType sum = storm::utility::zero(); + auto sum = storm::utility::zero(); for (size_t j = 0; j < lambdas.size(); ++j) { if (lambdas[j] != storm::utility::zero()) { sum += lambdas[j] * result_backup.at( @@ -117,12 +138,19 @@ namespace storm { currentValue += iter->second * sum; } - // Update the selected actions + // Update the selected actions TODO make this nicer if ((min && cc.isLess(storm::utility::zero(), chosenValue - currentValue)) || (!min && cc.isLess(storm::utility::zero(), currentValue - chosenValue))) { chosenValue = currentValue; chosenActionIndex = action; + } else if ((min && cc.isEqual(storm::utility::zero(), + chosenValue - currentValue)) || + (!min && + cc.isEqual(storm::utility::zero(), + currentValue - chosenValue))) { + chosenValue = currentValue; + chosenActionIndex = action; } // TODO tie breaker? } @@ -148,8 +176,8 @@ namespace storm { STORM_PRINT("Grid approximation took " << iteration << " iterations" << std::endl); beliefGrid.push_back(initialBelief); - beliefIsKnown.push_back( - target_observations.find(initialBelief.observation) != target_observations.end()); + beliefIsTarget.push_back( + targetObservations.find(initialBelief.observation) != targetObservations.end()); std::pair>, std::vector> temp = computeSubSimplexAndLambdas( initialBelief.probabilities, gridResolution); @@ -167,92 +195,197 @@ namespace storm { // Now onto the under-approximation + std::set visitedBelieves; + std::deque believesToBeExpanded; + std::map beliefStateMap; + std::vector> transitions; + std::vector targetStates; + + uint64_t stateId = 0; + beliefStateMap[initialBelief.id] = stateId; + ++stateId; + + // Expand the believes TODO capsuling + visitedBelieves.insert(initialBelief.id); + believesToBeExpanded.push_back(initialBelief.id); + while (!believesToBeExpanded.empty()) { + auto currentBeliefId = believesToBeExpanded.front(); + std::map transitionsInState; + STORM_LOG_DEBUG("Exploring Belief " << beliefList[currentBeliefId].observation << "||" + << beliefList[currentBeliefId].probabilities); + if (beliefIsTarget[currentBeliefId]) { + // add a self-loop to target states and save them + transitionsInState[beliefStateMap[currentBeliefId]] = storm::utility::one(); + targetStates.push_back(beliefStateMap[currentBeliefId]); + } else { + if (chosenActions.find(currentBeliefId) == chosenActions.end()) { + // If the current Belief is not part of the grid, we have not computed the action to choose yet + chosenActions[currentBeliefId] = extractBestAction(pomdp, beliefList, beliefIsTarget, + targetObservations, + observationProbabilities, + nextBelieves, result, gridResolution, + currentBeliefId, beliefList.size(), min); + } + for (auto iter = observationProbabilities[currentBeliefId][chosenActions[currentBeliefId]].begin(); + iter != + observationProbabilities[currentBeliefId][chosenActions[currentBeliefId]].end(); ++iter) { + uint32_t observation = iter->first; + uint64_t nextBeliefId = nextBelieves[currentBeliefId][chosenActions[currentBeliefId]][observation]; + if (visitedBelieves.insert(nextBeliefId).second) { + beliefStateMap[nextBeliefId] = stateId; + ++stateId; + believesToBeExpanded.push_back(nextBeliefId); + } + transitionsInState[beliefStateMap[nextBeliefId]] = iter->second; + } + } + transitions.push_back(transitionsInState); + believesToBeExpanded.pop_front(); + } + + for (size_t i = 0; i < transitions.size(); ++i) { + for (auto const &transition : transitions[i]) { + STORM_LOG_DEBUG( + "Transition: " << i << " -- " << transition.second << "--> " << transition.first); + } + } + storm::models::sparse::StateLabeling labeling(transitions.size()); + labeling.addLabel("init"); + labeling.addLabel("target"); + labeling.addLabelToState("init", 0); + for (auto targetState : targetStates) { + labeling.addLabelToState("target", targetState); + } + storm::storage::sparse::ModelComponents modelComponents( + buildTransitionMatrix(transitions), labeling); + + storm::models::sparse::Dtmc underApproxDtmc(modelComponents); + auto model = std::make_shared>(underApproxDtmc); + model->printModelInformationToStream(std::cout); + + std::string propertyString = min ? "Pmin=? [F \"target\"]" : "Pmax=? [F \"target\"]"; + std::vector propertyVector = storm::api::parseProperties(propertyString); + std::shared_ptr property = storm::api::extractFormulasFromProperties( + propertyVector).front(); + + std::unique_ptr res( + storm::api::verifyWithSparseEngine(model, storm::api::createTask(property, + true))); + STORM_LOG_ASSERT(res, "Result does not exist."); + res->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates())); + ValueType resultValue = res->asExplicitQuantitativeCheckResult().getValueMap().begin()->second; + + + STORM_PRINT("Over-Approximation Result: " << overApprox << std::endl); + STORM_PRINT("Under-Approximation Result: " << resultValue << std::endl); + } + + template + storm::storage::SparseMatrix + ApproximatePOMDPModelchecker::buildTransitionMatrix( + std::vector> transitions) { + uint_fast64_t currentRow = 0; + uint64_t nrEntries = 0; + for (auto const &map : transitions) { + nrEntries += map.size(); + } + storm::storage::SparseMatrixBuilder smb(transitions.size(), transitions.size(), nrEntries); + for (auto const &map : transitions) { + for (auto const &transition : map) { + smb.addNextValue(currentRow, transition.first, transition.second); + } + ++currentRow; + } + return smb.build(); + } - /* All this has to be put into a separate function as we have to repeat it for other believes not in the grid - // first compute some things for the initial belief + template + uint64_t ApproximatePOMDPModelchecker::extractBestAction( + storm::models::sparse::Pomdp const &pomdp, + std::vector> &beliefList, + std::vector &beliefIsTarget, + std::set &targetObservations, + std::map>> &observationProbabilities, + std::map>> &nextBelieves, + std::map &result, + uint64_t gridResolution, uint64_t currentBeliefId, uint64_t nextId, bool min) { + auto cc = storm::utility::ConstantsComparator(); + storm::pomdp::Belief currentBelief = beliefList[currentBeliefId]; std::vector> observationProbabilitiesInAction; std::vector> nextBelievesInAction; - uint64_t initialNumChoices = pomdp.getNumberOfChoices( - pomdp.getStatesWithObservation(initialBelief.observation).front()); - for (uint64_t action = 0; action < initialNumChoices; ++action) { + uint64_t numChoices = pomdp.getNumberOfChoices( + pomdp.getStatesWithObservation(currentBelief.observation).front()); + for (uint64_t action = 0; action < numChoices; ++action) { std::map actionObservationProbabilities = computeObservationProbabilitiesAfterAction( - pomdp, initialBelief, action); + pomdp, currentBelief, action); std::map actionObservationBelieves; for (auto iter = actionObservationProbabilities.begin(); iter != actionObservationProbabilities.end(); ++iter) { uint32_t observation = iter->first; actionObservationBelieves[observation] = getBeliefAfterActionAndObservation(pomdp, beliefList, - initialBelief, + beliefIsTarget, + targetObservations, + currentBelief, action, observation, nextId); + nextId = beliefList.size(); } observationProbabilitiesInAction.push_back(actionObservationProbabilities); nextBelievesInAction.push_back(actionObservationBelieves); } - observationProbabilities.push_back(observationProbabilitiesInAction); - nextBelieves.push_back(nextBelievesInAction); + //STORM_LOG_DEBUG("ID " << currentBeliefId << " add " << observationProbabilitiesInAction); + observationProbabilities.emplace(std::make_pair(currentBeliefId, observationProbabilitiesInAction)); + nextBelieves.emplace(std::make_pair(currentBeliefId, nextBelievesInAction)); - // do one step here to get the best action in the initial state + // choose the action which results in the value computed by the over-approximation ValueType chosenValue = min ? storm::utility::infinity() : -storm::utility::infinity(); uint64_t chosenActionIndex = std::numeric_limits::infinity(); ValueType currentValue; - for (uint64_t action = 0; action < initialNumChoices; ++action) { + for (uint64_t action = 0; action < numChoices; ++action) { currentValue = storm::utility::zero(); // simply change this for rewards? - for (auto iter = observationProbabilities[initialBelief.id][action].begin(); - iter != observationProbabilities[initialBelief.id][action].end(); ++iter) { + for (auto iter = observationProbabilities[currentBelief.id][action].begin(); + iter != observationProbabilities[currentBelief.id][action].end(); ++iter) { uint32_t observation = iter->first; - storm::pomdp::Belief nextBelief = beliefList[nextBelieves[initialBelief.id][action][observation]]; + storm::pomdp::Belief nextBelief = beliefList[nextBelieves[currentBelief.id][action][observation]]; // compute subsimplex and lambdas according to the Lovejoy paper to approximate the next belief - temp = computeSubSimplexAndLambdas( + auto temp = computeSubSimplexAndLambdas( nextBelief.probabilities, gridResolution); std::vector> subSimplex = temp.first; std::vector lambdas = temp.second; - ValueType sum = storm::utility::zero(); + auto sum = storm::utility::zero(); for (size_t j = 0; j < lambdas.size(); ++j) { if (lambdas[j] != storm::utility::zero()) { sum += lambdas[j] * result.at( - getBeliefIdInVector(beliefGrid, observation, subSimplex[j])); + getBeliefIdInVector(beliefList, observation, subSimplex[j])); } } currentValue += iter->second * sum; } - // Update the selected actions + // Update the selected actions TODO make this nicer if ((min && cc.isLess(storm::utility::zero(), chosenValue - currentValue)) || (!min && cc.isLess(storm::utility::zero(), currentValue - chosenValue))) { chosenValue = currentValue; chosenActionIndex = action; + } else if ((min && cc.isEqual(storm::utility::zero(), chosenValue - currentValue)) || + (!min && + cc.isEqual(storm::utility::zero(), currentValue - chosenValue))) { + chosenValue = currentValue; + chosenActionIndex = action; } } - chosenActions[initialBelief.id] = chosenActionIndex;*/ - - std::set exploredBelieves; - std::deque believesToBeExplored; - - exploredBelieves.insert(initialBelief.id); - believesToBeExplored.push_back(initialBelief.id); - while (!believesToBeExplored.empty()) { - auto currentBeliefId = believesToBeExplored.front(); - if (chosenActions.find(currentBeliefId) != chosenActions.end()) { - - } else { - - } - believesToBeExplored.pop_front(); - } - - - STORM_PRINT("Over-Approximation Result: " << overApprox << std::endl); + return chosenActionIndex; } + template uint64_t ApproximatePOMDPModelchecker::getBeliefIdInVector( std::vector> &grid, uint32_t observation, @@ -260,12 +393,10 @@ namespace storm { for (auto const &belief : grid) { if (belief.observation == observation) { if (belief.probabilities == probabilities) { - STORM_LOG_DEBUG("Found belief with id " << std::to_string(belief.id)); return belief.id; } } } - STORM_LOG_DEBUG("Did not find the belief in the grid"); return -1; } @@ -482,7 +613,8 @@ namespace storm { template uint64_t ApproximatePOMDPModelchecker::getBeliefAfterActionAndObservation( storm::models::sparse::Pomdp const &pomdp, - std::vector> &beliefList, storm::pomdp::Belief belief, + std::vector> &beliefList, std::vector &beliefIsTarget, + std::set &targetObservations, storm::pomdp::Belief belief, uint64_t actionIndex, uint32_t observation, uint64_t id) { std::vector distributionAfter(pomdp.getNumberOfStates(), storm::utility::zero()); for (uint64_t state = 0; state < pomdp.getNumberOfStates(); ++state) { @@ -508,6 +640,7 @@ namespace storm { return getBeliefIdInVector(beliefList, observation, distributionAfter); } else { beliefList.push_back(storm::pomdp::Belief{id, observation, distributionAfter}); + beliefIsTarget.push_back(targetObservations.find(observation) != targetObservations.end()); return id; } } diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h index e5636d5a9..46b5ee1e4 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.h @@ -1,22 +1,24 @@ #include -#include "storm/modelchecker/CheckTask.h" +#include "storm/api/storm.h" #include "storm/models/sparse/Pomdp.h" #include "storm/utility/logging.h" #include "storm-pomdp/storage/Belief.h" +#include "storm/storage/jani/Property.h" + namespace storm { namespace pomdp { namespace modelchecker { class POMDPCheckResult; - template> + template> class ApproximatePOMDPModelchecker { public: explicit ApproximatePOMDPModelchecker(); /*std::unique_ptr*/ void computeReachabilityProbability(storm::models::sparse::Pomdp const &pomdp, - std::set target_observations, bool min, + std::set targetObservations, bool min, uint64_t gridResolution); std::unique_ptr @@ -24,6 +26,29 @@ namespace storm { std::set target_observations, uint64_t gridResolution); private: + /** + * TODO + * @param pomdp + * @param beliefList + * @param observationProbabilities + * @param nextBelieves + * @param result + * @param gridResolution + * @param currentBeliefId + * @param nextId + * @param min + * @return + */ + uint64_t extractBestAction(storm::models::sparse::Pomdp const &pomdp, + std::vector> &beliefList, + std::vector &beliefIsTarget, + std::set &target_observations, + std::map>> &observationProbabilities, + std::map>> &nextBelieves, + std::map &result, + uint64_t gridResolution, uint64_t currentBeliefId, uint64_t nextId, + bool min); + /** * * @param pomdp @@ -82,10 +107,13 @@ namespace storm { * @return the resulting belief (observation and distribution) */ uint64_t - getBeliefAfterActionAndObservation(const models::sparse::Pomdp &pomdp, - std::vector> &beliefList, - storm::pomdp::Belief belief, - uint64_t actionIndex, uint32_t observation, uint64_t id); + getBeliefAfterActionAndObservation( + storm::models::sparse::Pomdp const &pomdp, + std::vector> &beliefList, + std::vector &beliefIsTarget, + std::set &targetObservations, + storm::pomdp::Belief belief, + uint64_t actionIndex, uint32_t observation, uint64_t id); /** * Helper method to get the next belief that results from a belief by performing an action @@ -108,6 +136,9 @@ namespace storm { */ uint64_t getBeliefIdInVector(std::vector> &grid, uint32_t observation, std::vector probabilities); + + storm::storage::SparseMatrix + buildTransitionMatrix(std::vector> transitions); }; }