From fa10087fba05b53c0a1417340c884fd0e359fffa Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 5 May 2020 13:00:26 +0200 Subject: [PATCH] Implemented triangulation in a more dynamic way. --- .../ApproximatePOMDPModelCheckerOptions.h | 1 + .../ApproximatePOMDPModelchecker.cpp | 8 +- src/storm-pomdp/storage/BeliefManager.h | 161 ++++++++++++------ 3 files changed, 112 insertions(+), 58 deletions(-) diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelCheckerOptions.h b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelCheckerOptions.h index 6d977a902..c24248467 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelCheckerOptions.h +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelCheckerOptions.h @@ -38,6 +38,7 @@ namespace storm { ValueType obsThresholdIncrementFactor = storm::utility::convertNumber(0.1); ValueType numericPrecision = storm::NumberTraits::IsExact ? storm::utility::zero() : storm::utility::convertNumber(1e-9); /// Used to decide whether two beliefs are equal + bool dynamicTriangulation = true; // Sets whether the triangulation is done in a dynamic way (yielding more precise triangulations) }; } } diff --git a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp index c8bc31f10..d72d56a6d 100644 --- a/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp +++ b/src/storm-pomdp/modelchecker/ApproximatePOMDPModelchecker.cpp @@ -185,7 +185,7 @@ namespace storm { if (options.discretize) { std::vector observationResolutionVector(pomdp.getNrObservations(), options.resolutionInit); - auto manager = std::make_shared(pomdp, options.numericPrecision); + auto manager = std::make_shared(pomdp, options.numericPrecision, options.dynamicTriangulation ? BeliefManagerType::TriangulationMode::Dynamic : BeliefManagerType::TriangulationMode::Static); if (rewardModelName) { manager->setRewardModel(rewardModelName); } @@ -205,7 +205,7 @@ namespace storm { } } if (options.unfold) { // Underapproximation (uses a fresh Belief manager) - auto manager = std::make_shared(pomdp, options.numericPrecision); + auto manager = std::make_shared(pomdp, options.numericPrecision, options.dynamicTriangulation ? BeliefManagerType::TriangulationMode::Dynamic : BeliefManagerType::TriangulationMode::Static); if (rewardModelName) { manager->setRewardModel(rewardModelName); } @@ -243,7 +243,7 @@ namespace storm { HeuristicParameters overApproxHeuristicPar; if (options.discretize) { // Setup and build first OverApproximation observationResolutionVector = std::vector(pomdp.getNrObservations(), options.resolutionInit); - overApproxBeliefManager = std::make_shared(pomdp, options.numericPrecision); + overApproxBeliefManager = std::make_shared(pomdp, options.numericPrecision, options.dynamicTriangulation ? BeliefManagerType::TriangulationMode::Dynamic : BeliefManagerType::TriangulationMode::Static); if (rewardModelName) { overApproxBeliefManager->setRewardModel(rewardModelName); } @@ -267,7 +267,7 @@ namespace storm { std::shared_ptr underApproximation; HeuristicParameters underApproxHeuristicPar; if (options.unfold) { // Setup and build first UnderApproximation - underApproxBeliefManager = std::make_shared(pomdp, options.numericPrecision); + underApproxBeliefManager = std::make_shared(pomdp, options.numericPrecision, options.dynamicTriangulation ? BeliefManagerType::TriangulationMode::Dynamic : BeliefManagerType::TriangulationMode::Static); if (rewardModelName) { underApproxBeliefManager->setRewardModel(rewardModelName); } diff --git a/src/storm-pomdp/storage/BeliefManager.h b/src/storm-pomdp/storage/BeliefManager.h index bc3f559d1..5a135a151 100644 --- a/src/storm-pomdp/storage/BeliefManager.h +++ b/src/storm-pomdp/storage/BeliefManager.h @@ -20,7 +20,12 @@ namespace storm { typedef boost::container::flat_set BeliefSupportType; typedef uint64_t BeliefId; - BeliefManager(PomdpType const& pomdp, BeliefValueType const& precision) : pomdp(pomdp), cc(precision, false) { + enum class TriangulationMode { + Static, + Dynamic + }; + + BeliefManager(PomdpType const& pomdp, BeliefValueType const& precision, TriangulationMode const& triangulationMode) : pomdp(pomdp), cc(precision, false), triangulationMode(triangulationMode) { beliefToIdMap.resize(pomdp.getNrObservations()); initialBeliefId = computeInitialBelief(); } @@ -288,65 +293,112 @@ namespace storm { } }; - Triangulation triangulateBelief(BeliefType belief, uint64_t resolution) { - STORM_LOG_ASSERT(assertBelief(belief), "Input belief for triangulation is not valid."); + void triangulateBeliefFreudenthal(BeliefType const& belief, uint64_t const& resolution, Triangulation& result) { + STORM_LOG_ASSERT(resolution != 0, "Invalid resolution: 0"); StateType numEntries = belief.size(); + auto convResolution = storm::utility::convertNumber(resolution); + // This is the Freudenthal Triangulation as described in Lovejoy (a whole lotta math) + // Probabilities will be triangulated to values in 0/N, 1/N, 2/N, ..., N/N + // Variable names are mostly based on the paper + // However, we speed this up a little by exploiting that belief states usually have sparse support (i.e. numEntries is much smaller than pomdp.getNumberOfStates()). + // Initialize diffs and the first row of the 'qs' matrix (aka v) + std::set> sorted_diffs; // d (and p?) in the paper + std::vector qsRow; // Row of the 'qs' matrix from the paper (initially corresponds to v + qsRow.reserve(numEntries); + std::vector toOriginalIndicesMap; // Maps 'local' indices to the original pomdp state indices + toOriginalIndicesMap.reserve(numEntries); + BeliefValueType x = convResolution; + for (auto const& entry : belief) { + qsRow.push_back(storm::utility::floor(x)); // v + sorted_diffs.emplace(toOriginalIndicesMap.size(), x - qsRow.back()); // x-v + toOriginalIndicesMap.push_back(entry.first); + x -= entry.second * convResolution; + } + // Insert a dummy 0 column in the qs matrix so the loops below are a bit simpler + qsRow.push_back(storm::utility::zero()); + + result.weights.reserve(numEntries); + result.gridPoints.reserve(numEntries); + auto currentSortedDiff = sorted_diffs.begin(); + auto previousSortedDiff = sorted_diffs.end(); + --previousSortedDiff; + for (StateType i = 0; i < numEntries; ++i) { + // Compute the weight for the grid points + BeliefValueType weight = previousSortedDiff->diff - currentSortedDiff->diff; + if (i == 0) { + // The first weight is a bit different + weight += storm::utility::one(); + } else { + // 'compute' the next row of the qs matrix + qsRow[previousSortedDiff->dimension] += storm::utility::one(); + } + if (!cc.isZero(weight)) { + result.weights.push_back(weight); + // Compute the grid point + BeliefType gridPoint; + for (StateType j = 0; j < numEntries; ++j) { + BeliefValueType gridPointEntry = qsRow[j] - qsRow[j + 1]; + if (!cc.isZero(gridPointEntry)) { + gridPoint[toOriginalIndicesMap[j]] = gridPointEntry / convResolution; + } + } + result.gridPoints.push_back(getOrAddBeliefId(gridPoint)); + } + previousSortedDiff = currentSortedDiff++; + } + } + + void triangulateBeliefDynamic(BeliefType const& belief, uint64_t const& resolution, Triangulation& result) { + // Find the best resolution for this belief, i.e., N such that the largest distance between one of the belief values to a value in {i/N | 0 ≤ i ≤ N} is minimal + uint64_t finalResolution = resolution; + BeliefValueType finalResolutionDist = storm::utility::one(); + // We don't need to check resolutions that are smaller than the maximal resolution divided by 2 (as we already checked multiples of these) + for (uint64_t currResolution = resolution; currResolution > resolution / 2; --currResolution) { + BeliefValueType currResDist = storm::utility::zero(); + BeliefValueType currResolutionConverted = storm::utility::convertNumber(currResolution); + bool continueWithNextResolution = false; + for (auto const& belEntry : belief) { + BeliefValueType product = belEntry.second * currResolutionConverted; + BeliefValueType dist = storm::utility::abs(product - storm::utility::round(product)) / currResolutionConverted; + if (dist > currResDist) { + if (dist > finalResolutionDist) { + // This resolution is worse than a previous resolution + continueWithNextResolution = true; + break; + } + currResDist = dist; + } + } + STORM_LOG_ASSERT(continueWithNextResolution || currResDist <= finalResolutionDist, "Distance for this resolution should not be larger than a previously checked one."); + if (!continueWithNextResolution) { + finalResolution = currResolution; + finalResolutionDist = currResDist; + } + } + + STORM_LOG_TRACE("Picking resolution " << finalResolution << " for belief " << toString(belief)); + + // do standard freudenthal with the found resolution + triangulateBeliefFreudenthal(belief, finalResolution, result); + } + + Triangulation triangulateBelief(BeliefType const& belief, uint64_t const& resolution) { + STORM_LOG_ASSERT(assertBelief(belief), "Input belief for triangulation is not valid."); Triangulation result; - // Quickly triangulate Dirac beliefs - if (numEntries == 1u) { + if (belief.size() == 1u) { result.weights.push_back(storm::utility::one()); result.gridPoints.push_back(getOrAddBeliefId(belief)); } else { - - auto convResolution = storm::utility::convertNumber(resolution); - // This is the Freudenthal Triangulation as described in Lovejoy (a whole lotta math) - // Variable names are mostly based on the paper - // However, we speed this up a little by exploiting that belief states usually have sparse support (i.e. numEntries is much smaller than pomdp.getNumberOfStates()). - // Initialize diffs and the first row of the 'qs' matrix (aka v) - std::set> sorted_diffs; // d (and p?) in the paper - std::vector qsRow; // Row of the 'qs' matrix from the paper (initially corresponds to v - qsRow.reserve(numEntries); - std::vector toOriginalIndicesMap; // Maps 'local' indices to the original pomdp state indices - toOriginalIndicesMap.reserve(numEntries); - BeliefValueType x = convResolution; - for (auto const& entry : belief) { - qsRow.push_back(storm::utility::floor(x)); // v - sorted_diffs.emplace(toOriginalIndicesMap.size(), x - qsRow.back()); // x-v - toOriginalIndicesMap.push_back(entry.first); - x -= entry.second * convResolution; - } - // Insert a dummy 0 column in the qs matrix so the loops below are a bit simpler - qsRow.push_back(storm::utility::zero()); - - result.weights.reserve(numEntries); - result.gridPoints.reserve(numEntries); - auto currentSortedDiff = sorted_diffs.begin(); - auto previousSortedDiff = sorted_diffs.end(); - --previousSortedDiff; - for (StateType i = 0; i < numEntries; ++i) { - // Compute the weight for the grid points - BeliefValueType weight = previousSortedDiff->diff - currentSortedDiff->diff; - if (i == 0) { - // The first weight is a bit different - weight += storm::utility::one(); - } else { - // 'compute' the next row of the qs matrix - qsRow[previousSortedDiff->dimension] += storm::utility::one(); - } - if (!cc.isZero(weight)) { - result.weights.push_back(weight); - // Compute the grid point - BeliefType gridPoint; - for (StateType j = 0; j < numEntries; ++j) { - BeliefValueType gridPointEntry = qsRow[j] - qsRow[j + 1]; - if (!cc.isZero(gridPointEntry)) { - gridPoint[toOriginalIndicesMap[j]] = gridPointEntry / convResolution; - } - } - result.gridPoints.push_back(getOrAddBeliefId(gridPoint)); - } - previousSortedDiff = currentSortedDiff++; + switch (triangulationMode) { + case TriangulationMode::Static: + triangulateBeliefFreudenthal(belief, resolution, result); + break; + case TriangulationMode::Dynamic: + triangulateBeliefDynamic(belief, resolution, result); + break; + default: + STORM_LOG_ASSERT(false, "Invalid triangulation mode."); } } STORM_LOG_ASSERT(assertTriangulation(belief, result), "Incorrect triangulation: " << toString(result)); @@ -445,6 +497,7 @@ namespace storm { storm::utility::ConstantsComparator cc; + TriangulationMode triangulationMode; }; }