implemented relative precision

src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp

@@ -135,7 +135,7 @@ namespace storm {
             }
 
             template <typename ModelType, typename GeometryValueType>
-            std::pair<std::vector<std::vector<GeometryValueType>>, std::vector<std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>>>> DeterministicSchedsLpChecker<ModelType, GeometryValueType>::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps) {
+            std::pair<std::vector<std::vector<GeometryValueType>>, std::vector<std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>>>> DeterministicSchedsLpChecker<ModelType, GeometryValueType>::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& 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<GeometryValueType>();
-                    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<ValueType>(milpGap), false);
                     gurobiLpModel->update();
                 }
@@ -644,7 +642,7 @@ namespace storm {
             }
             
             template <typename ModelType, typename GeometryValueType>
-            void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& infeasableAreas, uint64_t const& depth) {
+            void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, Point const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& 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<GeometryValueType>(gurobiLpModel->getContinuousValue(currentObjectiveVariables[objIndex], solutionIndex)));
-                                    if (p.back() > newPoint[objIndex] + eps / storm::utility::convertNumber<GeometryValueType>(2)) {
+                                    if (p.back() > newPoint[objIndex] + eps[objIndex] / storm::utility::convertNumber<GeometryValueType>(2)) {
                                         // The other solution dominates the newPoint in this dimension and is also not too close to the newPoint.
                                         considerP = true;
                                     }

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<Point>, std::vector<Polytope>> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps);
+                std::pair<std::vector<Point>, std::vector<Polytope>> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, Point const& eps);
 
             private:
                 std::vector<std::vector<storm::expressions::Expression>> 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<GeometryValueType>& polytopeTree, GeometryValueType const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& infeasableAreas, uint64_t const& depth);
+                void checkRecursive(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, Point const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& infeasableAreas, uint64_t const& depth);
                 
                 ModelType const& model;
                 std::vector<DeterministicSchedsObjectiveHelper<ModelType>> const& objectiveHelper;

src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp

@@ -294,6 +294,21 @@ namespace storm {
             std::unique_ptr<CheckResult> DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::check(Environment const& env) {
                 clean();
                 initializeFacets(env);
+                
+                // Compute the relative precision in each dimension.
+                // Todo: also allow for absolute precision
+                eps = std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>());
+                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<GeometryValueType>(env.modelchecker().multi().getPrecision());
+                epsScalingFactor += epsScalingFactor;
+                epsScalingFactor = epsScalingFactor / storm::utility::convertNumber<GeometryValueType, uint64_t>(pointset.size());
+                storm::utility::vector::scaleVectorInPlace(eps, epsScalingFactor);
+                std::cout << "scaled precision is " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(eps)) << std::endl;
+                
                 while (!unprocessedFacets.empty()) {
                     Facet f = std::move(unprocessedFacets.front());
                     unprocessedFacets.pop();
@@ -449,14 +464,11 @@ namespace storm {
 
             template <class SparseModelType, typename GeometryValueType>
             void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::processFacet(Environment const& env, Facet& f) {
-                GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(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 <class SparseModelType, typename GeometryValueType>
-            bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet& f, GeometryValueType const& eps) {
+            bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet& f) {
                 
                 // Invoke optimization and insert the explored points
                 boost::optional<PointId> optPointId;

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<GeometryValueType>& vector) const;
@@ -180,7 +180,7 @@ namespace storm {
                 std::queue<Facet> unprocessedFacets;
                 Polytope overApproximation;
                 std::vector<Polytope> unachievableAreas;
-                
+                std::vector<GeometryValueType> eps;
                 std::shared_ptr<DeterministicSchedsLpChecker<SparseModelType, GeometryValueType>> lpChecker;
                 std::unique_ptr<PcaaWeightVectorChecker<SparseModelType>> wvChecker;
                 std::vector<DeterministicSchedsObjectiveHelper<SparseModelType>> objectiveHelper;

src/storm/storage/geometry/PolytopeTree.h

@@ -70,20 +70,23 @@ 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<ValueType> const& point, ValueType const& offset = storm::utility::zero<ValueType>()) {
-                    if (storm::utility::isZero(offset)) {
+                void substractDownwardClosure(std::vector<ValueType> const& point) {
                     setMinus(Polytope<ValueType>::createDownwardClosure({point}));
-                    } else {
-                        std::vector<ValueType> pointPrime;
-                        pointPrime.reserve(point.size());
-                        for (auto const& coordinate : point) {
-                            pointPrime.push_back(coordinate + offset);
                 }
+                
+                /*!
+                 * 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<ValueType> const& point, std::vector<ValueType> const& offset) {
+                    assert(point.size() == offset.size());
+                    std::vector<ValueType> pointPrime(point.size());
+                    storm::utility::vector::addVectors(point, offset, pointPrime);
                     setMinus(Polytope<ValueType>::createDownwardClosure({pointPrime}));
                 }
-                }
                 
                 /*!
                  * Returns true if this is the empty set.