From 74f78102334160d002e78417d532ab4e3f4b4d8f Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Fri, 5 Jul 2019 16:57:57 +0200 Subject: [PATCH] implemented relative precision --- .../DeterministicSchedsLpChecker.cpp | 14 +++++----- .../DeterministicSchedsLpChecker.h | 4 +-- .../DeterministicSchedsParetoExplorer.cpp | 22 +++++++++++---- .../DeterministicSchedsParetoExplorer.h | 4 +-- src/storm/storage/geometry/PolytopeTree.h | 27 ++++++++++--------- 5 files changed, 42 insertions(+), 29 deletions(-) diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index 44aabf04e..6d261a829 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -135,7 +135,7 @@ namespace storm { } template - std::pair>, std::vector>>> DeterministicSchedsLpChecker::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps) { + std::pair>, std::vector>>> DeterministicSchedsLpChecker::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, Point const& eps) { STORM_LOG_INFO("Checking " << polytopeTree.toString()); swCheck.start(); STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); @@ -145,11 +145,9 @@ namespace storm { if (gurobiLpModel) { // For gurobi, it is possible to specify a gap between the obtained lower/upper objective bounds. - GeometryValueType milpGap = storm::utility::zero(); - for (auto const& wi : currentWeightVector) { - milpGap += wi; - } - milpGap *= eps; + // Let p be the found solution point, q be the optimal (unknown) solution point, and w be the current weight vector. + // The gap between the solution p and q is |w*p - w*q| = |w*(p-q)| + GeometryValueType milpGap = storm::utility::vector::dotProduct(currentWeightVector, eps); gurobiLpModel->setMaximalMILPGap(storm::utility::convertNumber(milpGap), false); gurobiLpModel->update(); } @@ -644,7 +642,7 @@ namespace storm { } template - void DeterministicSchedsLpChecker::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth) { + void DeterministicSchedsLpChecker::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, Point const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth) { std::cout << "."; std::cout.flush(); STORM_LOG_ASSERT(!polytopeTree.isEmpty(), "Tree node is empty"); @@ -693,7 +691,7 @@ namespace storm { bool considerP = false; for (uint64_t objIndex = 0; objIndex < currentObjectiveVariables.size(); ++objIndex) { p.push_back(storm::utility::convertNumber(gurobiLpModel->getContinuousValue(currentObjectiveVariables[objIndex], solutionIndex))); - if (p.back() > newPoint[objIndex] + eps / storm::utility::convertNumber(2)) { + if (p.back() > newPoint[objIndex] + eps[objIndex] / storm::utility::convertNumber(2)) { // The other solution dominates the newPoint in this dimension and is also not too close to the newPoint. considerP = true; } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h index aa04ad926..4416bb307 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h @@ -43,7 +43,7 @@ namespace storm { * Optimizes in the currently given direction, recursively checks for points in the given area. * @return all pareto optimal points in the area given by polytopeTree as well as a set of area in which no solution lies (the points might be achievable via some point outside of this area, though) */ - std::pair, std::vector> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps); + std::pair, std::vector> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, Point const& eps); private: std::vector> createEcVariables(); @@ -52,7 +52,7 @@ namespace storm { // Builds the induced markov chain of the current model and checks whether the resulting value coincide with the result of the lp solver. void validateCurrentModel(Environment const& env) const; - void checkRecursive(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth); + void checkRecursive(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, Point const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth); ModelType const& model; std::vector> const& objectiveHelper; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp index a1e05f377..21bc8fe7a 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp @@ -294,6 +294,21 @@ namespace storm { std::unique_ptr DeterministicSchedsParetoExplorer::check(Environment const& env) { clean(); initializeFacets(env); + + // Compute the relative precision in each dimension. + // Todo: also allow for absolute precision + eps = std::vector(objectives.size(), storm::utility::zero()); + for (auto const& point : pointset) { + for (uint64_t i = 0; i < eps.size(); ++i) { + eps[i] += storm::utility::abs(point.second.get()[i]); + } + } + GeometryValueType epsScalingFactor = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); + epsScalingFactor += epsScalingFactor; + epsScalingFactor = epsScalingFactor / storm::utility::convertNumber(pointset.size()); + storm::utility::vector::scaleVectorInPlace(eps, epsScalingFactor); + std::cout << "scaled precision is " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(eps)) << std::endl; + while (!unprocessedFacets.empty()) { Facet f = std::move(unprocessedFacets.front()); unprocessedFacets.pop(); @@ -449,14 +464,11 @@ namespace storm { template void DeterministicSchedsParetoExplorer::processFacet(Environment const& env, Facet& f) { - GeometryValueType eps = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); - eps += eps; // The unknown area (box) can actually have size 2*eps - if (!wvChecker) { lpChecker->setCurrentWeightVector(f.getHalfspace().normalVector()); } - if (optimizeAndSplitFacet(env, f, eps)) { + if (optimizeAndSplitFacet(env, f)) { return; } @@ -492,7 +504,7 @@ namespace storm { } template - bool DeterministicSchedsParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet& f, GeometryValueType const& eps) { + bool DeterministicSchedsParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet& f) { // Invoke optimization and insert the explored points boost::optional optPointId; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h index 7751d9bf8..68ea2511e 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h @@ -171,7 +171,7 @@ namespace storm { * 2. New facets are generated and (if not already precise enough) added to unprocessedFacets * 3. true is returned */ - bool optimizeAndSplitFacet(Environment const& env, Facet& f, GeometryValueType const& eps); + bool optimizeAndSplitFacet(Environment const& env, Facet& f); Polytope negateMinObjectives(Polytope const& polytope) const; void negateMinObjectives(std::vector& vector) const; @@ -180,7 +180,7 @@ namespace storm { std::queue unprocessedFacets; Polytope overApproximation; std::vector unachievableAreas; - + std::vector eps; std::shared_ptr> lpChecker; std::unique_ptr> wvChecker; std::vector> objectiveHelper; diff --git a/src/storm/storage/geometry/PolytopeTree.h b/src/storm/storage/geometry/PolytopeTree.h index 34a5b16a6..016c8cbcc 100644 --- a/src/storm/storage/geometry/PolytopeTree.h +++ b/src/storm/storage/geometry/PolytopeTree.h @@ -70,19 +70,22 @@ namespace storm { /*! * Substracts the downward closure of the given point from this set. * @param point the given point - * @param offset + * @param offset coordinates that are added to the point before taking its downward closure */ - void substractDownwardClosure(std::vector const& point, ValueType const& offset = storm::utility::zero()) { - if (storm::utility::isZero(offset)) { - setMinus(Polytope::createDownwardClosure({point})); - } else { - std::vector pointPrime; - pointPrime.reserve(point.size()); - for (auto const& coordinate : point) { - pointPrime.push_back(coordinate + offset); - } - setMinus(Polytope::createDownwardClosure({pointPrime})); - } + void substractDownwardClosure(std::vector const& point) { + setMinus(Polytope::createDownwardClosure({point})); + } + + /*! + * Substracts the downward closure of the given point from this set. + * @param point the given point + * @param offset coordinates that are added to the point before taking its downward closure + */ + void substractDownwardClosure(std::vector const& point, std::vector const& offset) { + assert(point.size() == offset.size()); + std::vector pointPrime(point.size()); + storm::utility::vector::addVectors(point, offset, pointPrime); + setMinus(Polytope::createDownwardClosure({pointPrime})); } /*!