From f1a57494e8eef741eb238b5acb2e71afa68d906e Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Fri, 3 May 2019 15:53:52 +0200
Subject: [PATCH] DeterministicSchedsParetoExplorer: Stop splitting facets if
 they are already 'small' enough.

---
 .../DeterministicSchedsParetoExplorer.cpp     | 34 +++++++++++++++----
 .../DeterministicSchedsParetoExplorer.h       |  2 +-
 2 files changed, 29 insertions(+), 7 deletions(-)

diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp
index 7306b34a6..09a7506e0 100644
--- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp
+++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp
@@ -10,6 +10,7 @@
 #include "storm/modelchecker/multiobjective/MultiObjectivePostprocessing.h"
 #include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h"
 #include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h"
+#include "storm/environment/solver/MinMaxSolverEnvironment.h"
 
 #include "storm/utility/export.h"
 #include "storm/utility/solver.h"
@@ -441,12 +442,13 @@ namespace storm {
 
             template <class SparseModelType, typename GeometryValueType>
             void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::processFacet(Environment const& env, Facet& f) {
-                if (optimizeAndSplitFacet(env,f)) {
+                GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
+                eps += eps; // The unknown area (box) can actually have size 2*eps
+
+                if (optimizeAndSplitFacet(env, f, eps)) {
                     return;
                 }
                 
-                GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
-                eps += eps; // The unknown area (box) can actually have size 2*eps
                 storm::storage::geometry::PolytopeTree<GeometryValueType> polytopeTree(f.getInducedPolytope(pointset, getReferenceCoordinates(env)));
                 for (auto const& point : pointset) {
                     polytopeTree.substractDownwardClosure(point.second.get(), eps);
@@ -466,8 +468,18 @@ namespace storm {
                 }
             }
             
+            template <typename GeometryValueType>
+            bool closePoints(std::vector<GeometryValueType> const& first, std::vector<GeometryValueType> const& second, GeometryValueType const& maxDistance) {
+                for (uint64_t i = 0; i < first.size(); ++i) {
+                    if (storm::utility::abs<GeometryValueType>(first[i] - second[i]) > maxDistance) {
+                        return false;
+                    }
+                }
+                return true;
+            }
+            
             template <class SparseModelType, typename GeometryValueType>
-            bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet& f) {
+            bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet& f, GeometryValueType const& eps) {
                 
                 // Invoke optimization and insert the explored points
                 boost::optional<PointId> optPointId;
@@ -489,17 +501,27 @@ namespace storm {
                 // Potentially generate new facets
                 if (optPointId) {
                     auto const& optPoint = pointset.getPoint(*optPointId);
-                    // TODO: this check might suffer from numerical errors. Check how much this would hurt us.
                     if (f.getHalfspace().contains(optPoint.get())) {
                         // The point is contained in the halfspace which means that no more splitting is possible.
                         return false;
                     } else {
-                        // Found a new Pareto optimal point -> generate new facets
+                        // Collect the new point with its neighbors.
+                        // Also check whether the remamining area is already sufficiently small.
+                        storm::storage::geometry::PolytopeTree<GeometryValueType> remainingArea(overApproximation->intersection(f.getHalfspace().invert()));
                         std::vector<std::vector<GeometryValueType>> vertices;
                         vertices.push_back(optPoint.get());
+                        auto minmaxPrec = storm::utility::convertNumber<GeometryValueType>(env.solver().minMax().getPrecision());
+                        minmaxPrec += minmaxPrec;
                         for (auto const& pId : f.getPoints()) {
                             vertices.push_back(pointset.getPoint(pId).get());
+                            remainingArea.substractDownwardClosure(vertices.back(), eps);
+                            STORM_LOG_WARN_COND((std::is_same<ModelValueType, storm::RationalNumber>::value) || !closePoints(vertices.front(), vertices.back(), minmaxPrec), "Found Pareto optimal points that are close to each other. This can be due to numerical issues. Maybe try exact mode?");
+                        }
+                        if (remainingArea.isEmpty()) {
+                            return false;
                         }
+                        
+                        // We need to generate new facets
                         auto newHalfspaceCandidates = storm::storage::geometry::Polytope<GeometryValueType>::createSelectiveDownwardClosure(vertices, storm::utility::vector::filterZero(f.getHalfspace().normalVector()))->getHalfspaces();
                         for (auto& h : newHalfspaceCandidates) {
                             if (!storm::utility::vector::hasNegativeEntry(h.normalVector())) {
diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h
index a3e3e0c1c..7751d9bf8 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);
+                bool optimizeAndSplitFacet(Environment const& env, Facet& f, GeometryValueType const& eps);
                 
                 Polytope negateMinObjectives(Polytope const& polytope) const;
                 void negateMinObjectives(std::vector<GeometryValueType>& vector) const;