implemented the LP solving to find a separating halfspace

Former-commit-id: d88558db0b

src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp

@@ -93,7 +93,46 @@ namespace storm {
                         return h;
                     }
                 }
-                return underApproximation->findSeparatingHalfspace(pointToBeSeparated);
+                
+                // We now use LP solving to find a seperating halfspace.
+                // Note that StoRM's LPSover does not support rational numbers. Hence, we optimize a polytope instead.
+                // The variables of the LP are the normalVector entries w_0 ... w_{n-1} and the minimal distance d between the halfspace and the paretoOptimalPoints
+                uint_fast64_t numVariables = pointToBeSeparated.size() + 1;
+                std::vector<storm::storage::geometry::Halfspace<RationalNumberType>> constraints;
+                constraints.reserve((numVariables-1) + 2 + paretoOptimalPoints.size());
+                // w_i >= 0 and d>= 0   <==>  -w_i <= 0 and -d <= 0
+                for(uint_fast64_t i = 0; i<numVariables; ++i) {
+                    WeightVector constraintCoefficients(numVariables, storm::utility::zero<RationalNumberType>());
+                    constraintCoefficients[i] = -storm::utility::one<RationalNumberType>();
+                    RationalNumberType constraintOffset = storm::utility::zero<RationalNumberType>();
+                    constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+                }
+                // sum_i w_i == 1   <==>   sum_i w_i <=1 and sum_i -w_i <= -1
+                {
+                    WeightVector constraintCoefficients(numVariables, storm::utility::one<RationalNumberType>());
+                    constraintCoefficients.back() = storm::utility::zero<RationalNumberType>();
+                    RationalNumberType constraintOffset = storm::utility::one<RationalNumberType>();
+                    constraints.emplace_back(constraintCoefficients, constraintOffset);
+                    storm::utility::vector::scaleVectorInPlace(constraintCoefficients, -storm::utility::one<RationalNumberType>());
+                    constraintOffset = -storm::utility::one<RationalNumberType>();
+                    constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+                }
+                // let q=pointToBeSeparated. For each paretoPoint x: q*w - x*w >= d  <==> (x-q) * w + d <= 0
+                for(auto const& paretoPoint : paretoOptimalPoints){
+                    WeightVector constraintCoefficients(numVariables-1);
+                    storm::utility::vector::subtractVectors(paretoPoint.first, pointToBeSeparated, constraintCoefficients);
+                    constraintCoefficients.push_back(storm::utility::one<RationalNumberType>());
+                    RationalNumberType constraintOffset = storm::utility::zero<RationalNumberType>();
+                    constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+                }
+                
+                WeightVector optimizationVector(numVariables);
+                optimizationVector.back() = storm::utility::one<RationalNumberType>(); // maximize d
+                std::pair<WeightVector, bool> optimizeResult = storm::storage::geometry::Polytope<RationalNumberType>::create(constraints)->optimize(optimizationVector);
+                STORM_LOG_THROW(optimizeResult.second, storm::exceptions::UnexpectedException, "There seems to be no seperating halfspace.");
+                optimizeResult.first.pop_back(); //drop the distance
+                RationalNumberType offset = storm::utility::vector::dotProduct(optimizeResult.first, pointToBeSeparated);
+                return storm::storage::geometry::Halfspace<RationalNumberType>(std::move(optimizeResult.first), std::move(offset));
             }
             
             template <class SparseModelType, typename RationalNumberType>

src/storage/geometry/HyproPolytope.cpp

@@ -148,7 +148,7 @@ namespace storm {
             std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::downwardClosure(boost::optional<Point> const& upperBounds) const {
                 if(this->isUniversal() || this->isEmpty()) {
                     // In these cases, the polytope does not change, i.e., we return a copy of this
-                    return std::make_shared<HyproPolytope<ValueType>>(internPolytope);
+                    return std::make_shared<HyproPolytope<ValueType>>(*this);
                 }
                 // Only keep the halfspaces where all entries of the normal vector are non-negative
                 std::vector<hypro::Halfspace<ValueType>> halfspaces;
@@ -184,6 +184,22 @@ namespace storm {
                 return std::make_shared<HyproPolytope<ValueType>>(HyproPolytopeType(std::move(halfspaces)));
             }
             
+            
+            template <typename ValueType>
+            std::pair<typename HyproPolytope<ValueType>::Point, bool> HyproPolytope<ValueType>::optimize(Point const& direction) const {
+                hypro::EvaluationResult<ValueType> evalRes = internPolytope.evaluate(storm::adapters::toHypro(direction));
+                switch (evalRes.errorCode) {
+                    case hypro::SOLUTION::FEAS:
+                        return std::make_pair(storm::adapters::fromHypro(evalRes.optimumValue), true);
+                    case hypro::SOLUTION::UNKNOWN:
+                        STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected eror code for Polytope evaluation");
+                        return std::make_pair(Point(), false);
+                    default:
+                        //solution is infinity or infeasible
+                        return std::make_pair(Point(), false);
+                }
+            }
+            
             template <typename ValueType>
             bool HyproPolytope<ValueType>::isHyproPolytope() const {
                 return true;

src/storage/geometry/HyproPolytope.h

@@ -107,6 +107,15 @@ namespace storm {
                  * @param upperBounds If given, this vector is considered for y (hence, max{x_i | x i \in P does not need to be computed)
                  */
                 virtual std::shared_ptr<Polytope<ValueType>> downwardClosure(boost::optional<Point> const& upperBounds = boost::none) const override;
+                
+                /*!
+                 * Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,
+                 * a point that maximizes dotPorduct(point, direction).
+                 * If such a point does not exist, the returned bool is false. There are two reasons for this:
+                 * - The polytope is empty
+                 * - The polytope is not bounded in the given direction
+                 */
+                virtual std::pair<Point, bool> optimize(Point const& direction) const override;
             
                 virtual bool isHyproPolytope() const override;
                 

src/storage/geometry/Polytope.cpp

@@ -129,9 +129,9 @@ namespace storm {
             }
             
             template <typename ValueType>
-            Halfspace<ValueType> Polytope<ValueType>::findSeparatingHalfspace(Point const& pointToBeSeparated) const {
+            std::pair<typename Polytope<ValueType>::Point, bool> Polytope<ValueType>::optimize(Point const& direction) const {
                 STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented.");
-                return Halfspace<ValueType>(Point(), storm::utility::zero<ValueType>());
+                return std::make_pair(Point(), false);
             }
             
             template <typename ValueType>

src/storage/geometry/Polytope.h

@@ -99,15 +99,13 @@ namespace storm {
                 virtual std::shared_ptr<Polytope<ValueType>> downwardClosure(boost::optional<Point> const& upperBounds = boost::none) const;
                 
                 /*!
-                 * Returns a halfspace h that separates this polytope from the given point p, i.e.,
-                 * - p lies on the border of h and
-                 * - for each x in this polytope, it holds that h contains x
-                 *
-                 * A halfspace with maximal distance to the polytope is preferred.
-                 *
-                 * @param pointToBeSeparated the point that is to be seperated
+                 * Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,
+                 * a point that maximizes dotPorduct(point, direction).
+                 * If such a point does not exist, the returned bool is false. There are two reasons for this:
+                 * - The polytope is empty
+                 * - The polytope is not bounded in the given direction
                  */
-                virtual Halfspace<ValueType> findSeparatingHalfspace(Point const& pointToBeSeparated) const;
+                virtual std::pair<Point, bool> optimize(Point const& direction) const;
                 
                 /*
                  * Returns a string representation of this polytope.