First version of MILP based deterministic scheduler technique.

src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsSimplexChecker.h

@@ -0,0 +1,387 @@
+#pragma once
+
+#include <vector>
+
+#include "storm/storage/geometry/Polytope.h"
+#include "storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h"
+#include "storm/storage/expressions/Expressions.h"
+#include "storm/utility/solver.h"
+#include "storm/solver/LpSolver.h"
+#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "storm/modelchecker/csl/SparseCtmcCslModelChecker.h"
+#include "storm/storage/BitVector.h"
+#include "storm/utility/graph.h"
+#include "storm/utility/vector.h"
+#include "storm/utility/Stopwatch.h"
+
+namespace storm {
+    
+    class Environment;
+    
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            /*!
+             * Represents a set of points in euclidean space.
+             * The set is defined as the union of the polytopes at the leafs of the tree.
+             * The polytope at inner nodes is always the convex union of its children.
+             * The sets described by the children of a node are disjoint.
+             * A child is always non-empty, i.e., isEmpty() should only hold for the root node.
+             * @tparam GeometryValueType
+             */
+            template <typename GeometryValueType>
+            class PolytopeTree {
+                typedef typename std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> Polytope;
+                typedef typename std::vector<GeometryValueType> Point;
+
+            public:
+                PolytopeTree(Polytope const& polytope) : polytope(polytope) {}
+                
+                /*!
+                 * Substracts the given rhs from this polytope.
+                 */
+                void setMinus(Polytope const& rhs) {
+                    // This operation only has an effect if the intersection of this and rhs is non-empty.
+                    if (!isEmpty() && !polytope->intersection(rhs)->isEmpty()) {
+                        if (children.empty()) {
+                            // This is a leaf node.
+                            // Apply splitting.
+                            auto newChildren = polytope->setMinus(rhs);
+                            if (newChildren.empty()) {
+                                // Delete this node.
+                                polytope = nullptr;
+                            } else if (newChildren.size() == 1) {
+                                // Replace this node with its only child
+                                polytope = newChildren.front()->clean();
+                            } else {
+                                // Add the new children to this node. There is no need to traverse them.
+                                for (auto& c : newChildren) {
+                                    children.push_back(c->clean());
+                                }
+                            }
+                        } else {
+                            // This is an inner node. Traverse the children and set this to the convex union of its children.
+                            std::vector<PolytopeTree<GeometryValueType>> newChildren;
+                            Polytope newPolytope = nullptr;
+                            for (auto& c : children) {
+                                c.setMinus(rhs);
+                                if (c.polytope != nullptr) {
+                                    newChildren.push_back(c);
+                                    if (newPolytope) {
+                                        newPolytope->convexUnion(c.polytope);
+                                    } else {
+                                        newPolytope = c.polytope;
+                                    }
+                                }
+                            }
+                            polytope = newPolytope; // nullptr, if no children left
+                            children = std::move(newChildren);
+                        }
+                    }
+                }
+                
+                void substractDownwardClosure(Point const& point, GeometryValueType const& eps) {
+                    std::vector<GeometryValueType>(pointPlusEps);
+                    for (auto const& coordinate : point) {
+                        pointPlusEps.push_back(coordinate + eps);
+                    }
+                    auto downwardOfPoint = storm::storage::geometry::Polytope<GeometryValueType>::createDownwardClosure({pointPlusEps});
+                    setMinus(downwardOfPoint);
+                }
+
+                bool isEmpty() const {
+                    return polytope == nullptr;
+                }
+                
+                void clear() {
+                    children.clear();
+                    polytope = nullptr;
+                }
+                
+                Polytope getPolytope() const {
+                    return polytope;
+                }
+                
+                std::vector<PolytopeTree>& getChildren() {
+                    return children;
+                }
+                
+                std::string toString() {
+                    std::stringstream s;
+                    s << "PolytopeTree node with " << getChildren().size() << " children: " << getPolytope()->toString(true) << std::endl << "Vertices: ";
+                    auto vertices = getPolytope()->getVertices();
+                    for (auto const& v : vertices) {
+                        s << "[";
+                        for (auto const& vi : v) {
+                            s << storm::utility::convertNumber<double>(vi) << ",";
+                        }
+                        s << "]\t";
+                    }
+                    s << std::endl;
+                    return s.str();
+                }
+                
+            private:
+                
+                
+                
+                Polytope polytope;
+                std::vector<PolytopeTree<GeometryValueType>> children;
+            };
+            
+            template <typename ModelType, typename GeometryValueType>
+            class DetSchedsSimplexChecker {
+            public:
+                
+                typedef typename ModelType::ValueType ValueType;
+                typedef typename std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> Polytope;
+                typedef typename std::vector<GeometryValueType> Point;
+                
+                DetSchedsSimplexChecker(std::shared_ptr<MultiObjectiveSchedulerEvaluator<ModelType>> const& schedulerEvaluator) : schedulerEvaluator(schedulerEvaluator) {
+                    init();
+                }
+                
+                ~DetSchedsSimplexChecker() {
+                    std::cout << "SIMPLEX CHECKER: " << swInit << " seconds for initialization" << std::endl;
+                    std::cout << "SIMPLEX CHECKER: " << swCheck << " seconds for checking, including" << std::endl;
+                    std::cout << "\t " << swLpBuild << " seconds for LP building" << std::endl;
+                    std::cout << "\t " << swLpSolve << " seconds for LP solving" << std::endl;
+                    std::cout << "SIMPLEX CHECKER: " << swAux << " seconds for aux stuff" << std::endl;
+                }
+                
+                std::pair<std::vector<Point>, std::vector<Polytope>> check(storm::Environment const& env, std::vector<GeometryValueType> const& weightVector, PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps) {
+                    std::cout << "Checking a Simplex with weight vector " << storm::utility::vector::toString(weightVector) << std::endl << " and root " << polytopeTree.toString() << std::endl << "\t";
+                    if (polytopeTree.isEmpty()) {
+                        return {{}, {}};
+                    }
+                    swCheck.start();
+                    
+                    swLpBuild.start();
+                    lpModel->push();
+                    currentObjectiveVariables.clear();
+                    
+                    // set up objective function for the given weight vector
+                    for (uint64_t objIndex = 0; objIndex < initialStateResults.size(); ++objIndex) {
+                        currentObjectiveVariables.push_back(lpModel->addUnboundedContinuousVariable("w_" + std::to_string(objIndex), storm::utility::convertNumber<ValueType>(weightVector[objIndex])));
+                        lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == initialStateResults[objIndex]);
+                    }
+                    lpModel->update();
+                    swLpBuild.stop();
+
+                    auto result = checkRecursive(weightVector, polytopeTree, eps);
+                    
+                    swLpBuild.start();
+                    lpModel->pop();
+                    lpModel->update();
+                    swLpBuild.stop();
+                    swCheck.stop();
+                    std::cout << " done!" << std::endl;
+                    return result;
+                }
+                
+            private:
+                
+                std::pair<std::vector<Point>, std::vector<Polytope>> checkRecursive(std::vector<GeometryValueType> const& weightVector, PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps) {
+                    std::cout << ".";
+                    std::cout.flush();
+                    STORM_LOG_ASSERT(!polytopeTree.isEmpty(), "Tree node is empty");
+                    STORM_LOG_ASSERT(!polytopeTree.getPolytope()->isEmpty(), "Tree node is empty.");
+                    STORM_LOG_TRACE("Checking " << polytopeTree.toString());
+                    
+                    auto vertices = polytopeTree.getPolytope()->getVertices();
+                    
+                    std::vector<Point> foundPoints;
+                    std::vector<Polytope> infeasableAreas;
+                    
+                    swLpBuild.start();
+                    lpModel->push();
+                    // Assert the constraints of the current polytope
+                    auto nodeConstraints = polytopeTree.getPolytope()->getConstraints(lpModel->getManager(), currentObjectiveVariables);
+                    for (auto const& constr : nodeConstraints) {
+                        lpModel->addConstraint("", constr);
+                    }
+                    lpModel->update();
+                    swLpBuild.stop();
+                    
+                    if (polytopeTree.getChildren().empty()) {
+                        // At leaf nodes we need to perform the actual check.
+                        swLpSolve.start();
+                        lpModel->optimize();
+                        swLpSolve.stop();
+                        
+                        if (lpModel->isInfeasible()) {
+                            infeasableAreas.push_back(polytopeTree.getPolytope());
+                            polytopeTree.clear();
+                        } else {
+                            STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded.");
+                            Point newPoint;
+                            for (auto const& objVar : currentObjectiveVariables) {
+                                newPoint.push_back(storm::utility::convertNumber<GeometryValueType>(lpModel->getContinuousValue(objVar)));
+                            }
+                            auto halfspace = storm::storage::geometry::Halfspace<GeometryValueType>(weightVector, storm::utility::vector::dotProduct(weightVector, newPoint)).invert();
+                            infeasableAreas.push_back(polytopeTree.getPolytope()->intersection(halfspace));
+                            if (infeasableAreas.back()->isEmpty()) {
+                                infeasableAreas.pop_back();
+                            }
+                            swAux.start();
+                            polytopeTree.setMinus(storm::storage::geometry::Polytope<GeometryValueType>::create({halfspace}));
+                            foundPoints.push_back(newPoint);
+                            polytopeTree.substractDownwardClosure(newPoint, eps);
+                            swAux.stop();
+                            if (!polytopeTree.isEmpty()) {
+                                auto childRes = checkRecursive(weightVector, polytopeTree, eps);
+                                foundPoints.insert(foundPoints.end(), childRes.first.begin(), childRes.first.end());
+                                infeasableAreas.insert(infeasableAreas.end(), childRes.second.begin(), childRes.second.end());
+                            }
+                        }
+                    } else {
+                        // Traverse all the children.
+                        for (uint64_t childId = 0; childId < polytopeTree.getChildren().size(); ++childId) {
+                            auto childRes = checkRecursive(weightVector, polytopeTree.getChildren()[childId], eps);
+                            STORM_LOG_ASSERT(polytopeTree.getChildren()[childId].isEmpty(), "expected empty children.");
+                            // Make the results known to the right siblings
+                            for (auto const& newPoint : childRes.first) {
+                                for (uint64_t siblingId = childId + 1; siblingId < polytopeTree.getChildren().size(); ++siblingId) {
+                                    polytopeTree.getChildren()[siblingId].substractDownwardClosure(newPoint, eps);
+                                }
+                            }
+                            foundPoints.insert(foundPoints.end(), childRes.first.begin(), childRes.first.end());
+                            infeasableAreas.insert(infeasableAreas.end(), childRes.second.begin(), childRes.second.end());
+                        }
+                        // All children are empty now, so this becomes empty.
+                        polytopeTree.clear();
+                    }
+                    swLpBuild.start();
+                    lpModel->pop();
+                    swLpBuild.stop();
+                    return {foundPoints, infeasableAreas};
+                }
+
+                /* Todo
+                ValueType getChoiceValueSummand(Objective<ValueType> const& objective, uint64_t choiceIndex) {
+                    auto const& model = schedulerEvaluator->getModel();
+                    storm::modelchecker::SparsePropositionalModelChecker<ModelType> mc(model);
+                    auto const& formula = *objective.formula;
+                    if (formula.isProbabilityOperatorFormula() && formula.getSubformula().isUntilFormula()) {
+                        return storm::utility::zero<ValueType>();
+                    } else if (formula.getSubformula().isEventuallyFormula() && (formula.isRewardOperatorFormula() || formula.isTimeOperatorFormula())) {
+                        
+                        storm::storage::BitVector rew0States = mc.check(formula.getSubformula().asEventuallyFormula().getSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector();
+                        if (formula.isRewardOperatorFormula()) {
+                            auto const& rewModel = formula.asRewardOperatorFormula().hasRewardModelName() ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()) : model.getUniqueRewardModel();
+                            storm::storage::BitVector statesWithoutReward = rewModel.getStatesWithZeroReward(model.getTransitionMatrix());
+                            rew0States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), statesWithoutReward, rew0States);
+                        }
+                        storm::utility::vector::setVectorValues(results[objIndex], rew0States, storm::utility::zero<ValueType>());
+                        schedulerIndependentStates.push_back(std::move(rew0States));
+                    } else if (formula.isRewardOperatorFormula() && formula.getSubformula().isTotalRewardFormula()) {
+                        auto const& rewModel = formula.asRewardOperatorFormula().hasRewardModelName() ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()) : model.getUniqueRewardModel();
+                        storm::storage::BitVector statesWithoutReward = rewModel.getStatesWithZeroReward(model.getTransitionMatrix());
+                        storm::storage::BitVector rew0States = storm::utility::graph::performProbGreater0E(model.getBackwardTransitions(), statesWithoutReward, ~statesWithoutReward);
+                        rew0States.complement();
+                        storm::utility::vector::setVectorValues(results[objIndex], rew0States, storm::utility::zero<ValueType>());
+                        schedulerIndependentStates.push_back(std::move(rew0States));
+                    } else {
+                        STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported.");
+                    }
+                }*/
+                
+                void init() {
+                    swInit.start();
+                    auto const& model = schedulerEvaluator->getModel();
+                    auto const& objectives = schedulerEvaluator->getObjectives();
+                    uint64_t numStates = model.getNumberOfStates();
+                    lpModel = storm::utility::solver::getLpSolver<ValueType>("model");
+                    lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize);
+                    initialStateResults.clear();
+                    
+                    auto one = lpModel->getConstant(storm::utility::one<ValueType>());
+                    
+                    // Create choice variables and assert that at least one choice is taken at each state.
+                    std::vector<storm::expressions::Expression> choiceVars;
+                    choiceVars.reserve(model.getNumberOfChoices());
+                    for (uint64_t state = 0; state < numStates; ++state) {
+                        uint64_t numChoices = model.getNumberOfChoices(state);
+                        if (numChoices == 1) {
+                            choiceVars.emplace_back();
+                        } else {
+                            std::vector<storm::expressions::Expression> localChoices;
+                            for (uint64_t choice = 0; choice < numChoices; ++choice) {
+                                localChoices.push_back(lpModel->addBoundedIntegerVariable("c" + std::to_string(state) + "_" + std::to_string(choice), 0, 1).getExpression());
+                            }
+                            lpModel->addConstraint("", storm::expressions::sum(localChoices).reduceNesting() >= one);
+                            choiceVars.insert(choiceVars.end(), localChoices.begin(), localChoices.end());
+                        }
+                    }
+                    
+                    for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                        Objective<ValueType> const& objective = objectives[objIndex];
+                        storm::storage::BitVector const& schedulerIndependentStates = schedulerEvaluator->getSchedulerIndependentStates(objIndex);
+                        // Create state variables
+                        std::vector<storm::expressions::Expression> stateVars;
+                        stateVars.reserve(numStates);
+                        for (uint64_t state = 0; state < numStates; ++state) {
+                            if (schedulerIndependentStates.get(state)) {
+                                stateVars.push_back(lpModel->getConstant(schedulerEvaluator->getSchedulerIndependentStateResult(objIndex, state)));
+                            } else {
+                                stateVars.push_back(lpModel->addContinuousVariable("x" + std::to_string(objIndex) + "_" + std::to_string(state), objective.lowerResultBound, objective.upperResultBound).getExpression());
+                            }
+                            if (state == *model.getInitialStates().begin()) {
+                                initialStateResults.push_back(stateVars.back());
+                            }
+                        }
+                        
+                        // Create and assert choice values
+                        for (uint64_t state = 0; state < numStates; ++state) {
+                            if (schedulerIndependentStates.get(state)) {
+                                continue;
+                            }
+                            storm::expressions::Expression stateValue;
+                            uint64_t numChoices = model.getNumberOfChoices(state);
+                            for (uint64_t choice = 0; choice < numChoices; ++choice) {
+                                storm::expressions::Expression choiceValue;
+                                if (objective.formula)
+                                for (auto const& transition : model.getTransitionMatrix().getRow(state, choice)) {
+                                    storm::expressions::Expression transitionValue = lpModel->getConstant(transition.getValue()) * stateVars[transition.getColumn()];
+                                    if (choiceValue.isInitialized()) {
+                                        choiceValue = choiceValue + transitionValue;
+                                    } else {
+                                        choiceValue = transitionValue;
+                                    }
+                                }
+                                choiceValue = choiceValue.simplify().reduceNesting();
+                                if (numChoices == 1) {
+                                    lpModel->addConstraint("", stateVars[state] == choiceValue);
+                                } else {
+                                    uint64_t globalChoiceIndex = model.getTransitionMatrix().getRowGroupIndices()[state] + choice;
+                                    storm::expressions::Expression maxDiff = lpModel->getConstant(objective.upperResultBound.get()) * (one - choiceVars[globalChoiceIndex]);
+                                    lpModel->addConstraint("", stateVars[state] - choiceValue <= maxDiff);
+                                    lpModel->addConstraint("", choiceValue - stateVars[state] <= maxDiff);
+                                }
+                            }
+                        }
+                    }
+                    lpModel->update();
+                    swInit.stop();
+                }
+                
+                std::shared_ptr<MultiObjectiveSchedulerEvaluator<ModelType>> schedulerEvaluator;
+
+                std::unique_ptr<storm::solver::LpSolver<ValueType>> lpModel;
+                std::vector<storm::expressions::Expression> initialStateResults;
+                std::vector<storm::expressions::Variable> currentObjectiveVariables;
+
+                
+                storm::utility::Stopwatch swInit;
+                storm::utility::Stopwatch swCheck;
+                storm::utility::Stopwatch swLpSolve;
+                storm::utility::Stopwatch swLpBuild;
+                storm::utility::Stopwatch swAux;
+            };
+            
+        }
+    }
+}

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

@@ -54,6 +54,7 @@ namespace storm {
                                     choiceValue += weightVector[objIndex] * stateResults[objIndex][state];
                                 } else {
                                     ValueType objValue = storm::utility::zero<ValueType>();
+                                    // TODO: Choice rewards?
                                     for (auto const& entry : transitionMatrix.getRow(choiceOffset + choice)) {
                                         objValue += entry.getValue() * stateResults[objIndex][entry.getColumn()];
                                     }

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

@@ -146,7 +146,6 @@ namespace storm {
             
                 // Find dominated and dominating points
                 auto pointsIt = points.begin();
-                auto pointsItE = points.end();
                 while (pointsIt != points.end()) {
                     switch (point.getDominance(pointsIt->second)) {
                         case Point::DominanceResult::Incomparable:
@@ -280,25 +279,12 @@ namespace storm {
             
             
             template <class SparseModelType, typename GeometryValueType>
-            typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Polytope const& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::getInducedSimplex(Pointset const& pointset) {
+            typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Polytope const& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::getInducedSimplex(Pointset const& pointset, std::vector<GeometryValueType> const& referenceCoordinates) {
                 if (!inducedSimplex) {
-                    std::vector<std::vector<GeometryValueType>> vertices;
-                    STORM_LOG_ASSERT(getNumberOfPoints() > 0, "Tried to compute the induced simplex, but not enough points are given.");
-                    auto pointIdIt = paretoPointsOnFacet.begin();
-                    auto pointIdItE = paretoPointsOnFacet.end();
-                    vertices.push_back(pointset.getPoint(*pointIdIt).get());
-                    std::vector<GeometryValueType> minPoint = vertices.back();
-                    
-                    for (++pointIdIt; pointIdIt != pointIdItE; ++pointIdIt) {
-                        vertices.push_back(pointset.getPoint(*pointIdIt).get());
-                        auto pIt = vertices.back().begin();
-                        for (auto& mi : minPoint) {
-                            mi = std::min(mi, *pIt);
-                            ++pIt;
-                        }
+                    std::vector<std::vector<GeometryValueType>> vertices = {referenceCoordinates};
+                    for (auto const& pId : paretoPointsOnFacet) {
+                        vertices.push_back(pointset.getPoint(pId).get());
                     }
-                    vertices.push_back(std::move(minPoint));
-                    
                     // This facet might lie at the 'border', which means that the downward closure has to be taken in some directions
                     storm::storage::BitVector dimensionsForDownwardClosure = storm::utility::vector::filterZero(this->halfspace.normalVector());
                     STORM_LOG_ASSERT(dimensionsForDownwardClosure.getNumberOfSetBits() + vertices.size() >= halfspace.normalVector().size() + 1, "The number of points on the facet is insufficient");
@@ -321,6 +307,7 @@ namespace storm {
                 originalModelInitialState = *preprocessorResult.originalModel.getInitialStates().begin();
                 schedulerEvaluator = std::make_shared<MultiObjectiveSchedulerEvaluator<SparseModelType>>(preprocessorResult);
                 weightVectorChecker = std::make_shared<DetSchedsWeightVectorChecker<SparseModelType>>(schedulerEvaluator);
+                simplexChecker = std::make_shared<DetSchedsSimplexChecker<SparseModelType, GeometryValueType>>(schedulerEvaluator);
             }
 
             template <class SparseModelType, typename GeometryValueType>
@@ -417,7 +404,7 @@ namespace storm {
                     auto points = weightVectorChecker->check(env, weightVector);
                     bool last = true;
                     for (auto pIt = points.rbegin(); pIt != points.rend(); ++pIt) {
-                        for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                        for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                             if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
                                 (*pIt)[objIndex] *= -storm::utility::one<ModelValueType>();
                             }
@@ -450,8 +437,22 @@ namespace storm {
                 }
             }
             
+            template <class SparseModelType, typename GeometryValueType>
+            std::vector<GeometryValueType>  DeterministicParetoExplorer<SparseModelType, GeometryValueType>::getReferenceCoordinates() const {
+                std::vector<GeometryValueType> result;
+                for (auto const& obj : schedulerEvaluator->getObjectives()) {
+                    ModelValueType value = storm::solver::minimize(obj.formula->getOptimalityType()) ? obj.upperResultBound.get() : obj.lowerResultBound.get();
+                    result.push_back(storm::utility::convertNumber<GeometryValueType>(value));
+                }
+                return result;
+            }
+
+            
             template <class SparseModelType, typename GeometryValueType>
             bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::checkFacetPrecision(Environment const& env, Facet& f) {
+                // TODO:
+                return false;
+                /*
                 auto const& inducedSimplex = f.getInducedSimplex(pointset);
                 
                 GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
@@ -462,11 +463,12 @@ namespace storm {
                 // If the intersection of both polytopes is empty, it means that there can not be a point y in the simplex
                 // such that y-eps is also in the simplex, i.e., the facet is already precise enough.
                 return inducedSimplex->intersection(shiftedSimplex)->isEmpty();
+                 */
             }
             
             template <class SparseModelType, typename GeometryValueType>
             bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::checkFacetPrecision(Environment const& env, Facet& f, std::set<PointId> const& collectedSimplexPoints) {
-                // TODO
+                assert(false);
                 return false;
             }
             
@@ -476,6 +478,26 @@ namespace storm {
                     return;
                 }
                 
+                GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
+                eps += eps; // The unknown area (box) can actually have size 2*eps
+                PolytopeTree<GeometryValueType> polytopeTree(f.getInducedSimplex(pointset, getReferenceCoordinates()));
+                for (auto const& point : pointset) {
+                    polytopeTree.substractDownwardClosure(point.second.get(), eps);
+                    if (polytopeTree.isEmpty()) {
+                        break;
+                    }
+                }
+                if (!polytopeTree.isEmpty()) {
+                    auto res = simplexChecker->check(env, f.getHalfspace().normalVector(), polytopeTree, eps);
+                    for (auto const& infeasableArea : res.second) {
+                        addUnachievableArea(env, infeasableArea);
+                    }
+                    for (auto& achievablePoint : res.first) {
+                        pointset.addPoint(env, Point(std::move(achievablePoint)));
+                    }
+                }
+                
+                /*
                 FacetAnalysisContext context = createAnalysisContext(env, f);
                 
                 if (findAndCheckCachedPoints(env, context)) {
@@ -489,15 +511,16 @@ namespace storm {
                 if (analyzePointsInSimplex(env, context)) {
                     return;
                 }
-                
+                */
                 // Reaching this point means that the facet could not be analyzed completely.
-                STORM_LOG_ERROR("Facet " << f.getHalfspace().toString(true) << " could not be analyzed completely.");
+                //STORM_LOG_ERROR("Facet " << f.getHalfspace().toString(true) << " could not be analyzed completely.");
             }
             
             template <class SparseModelType, typename GeometryValueType>
             typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::FacetAnalysisContext DeterministicParetoExplorer<SparseModelType, GeometryValueType>::createAnalysisContext(Environment const& env, Facet& f) {
-                FacetAnalysisContext res(f);
                 
+                FacetAnalysisContext res(f);
+                /*
                 res.expressionManager = std::make_shared<storm::expressions::ExpressionManager>();
                 res.smtSolver = storm::utility::solver::SmtSolverFactory().create(*res.expressionManager);
                 
@@ -524,9 +547,7 @@ namespace storm {
                     }
                 }
                 res.smtSolver->add(xme);
-                
-                res.smtSolver->push();
-                
+                */
                 return res;
             }
             
@@ -535,7 +556,7 @@ namespace storm {
                 // Obtain the correct weight vector
                 auto weightVector = storm::utility::vector::convertNumericVector<ModelValueType>(f.getHalfspace().normalVector());
                 bool weightVectorYieldsParetoOptimalPoint = !storm::utility::vector::hasZeroEntry(weightVector);
-                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
                         weightVector[objIndex] *= -storm::utility::one<ModelValueType>();
                     }
@@ -546,7 +567,7 @@ namespace storm {
                 auto points = weightVectorChecker->check(env, weightVector);
                 bool last = true;
                 for (auto pIt = points.rbegin(); pIt != points.rend(); ++pIt) {
-                    for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                    for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                         if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
                             (*pIt)[objIndex] *= -storm::utility::one<ModelValueType>();
                         }
@@ -606,50 +627,80 @@ namespace storm {
             }
             
             template <class SparseModelType, typename GeometryValueType>
-            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::findAndCheckCachedPoints(Environment const& env, FacetAnalysisContext& context) {
-                Polytope inducedPoly = context.facet.getInducedSimplex(pointset);
-                pointset.collectPointsInPolytope(context.collectedPoints, inducedPoly);
-                if (context.collectedPoints.size() != context.facet.getNumberOfPoints()) {
-                    STORM_LOG_ASSERT(context.collectedPoints.size() > context.facet.getNumberOfPoints(), "Did not find all points on the facet");
-                    // formulate SMT constraints to assert that xMinusEps is not dominated by one of the cached points
-                    for (auto const& pId : context.collectedPoints) {
-                        auto const& coordinates = pointset.getPoint(pId).get();
-                        storm::expressions::Expression pointAchievesXMinusEps;
-                        for (uint64_t i = 0; i < coordinates.size(); ++i) {
-                            storm::expressions::Expression subExpr = context.xMinusEps[i] <= context.expressionManager->rational(coordinates[i]);
-                            if (i == 0) {
-                                pointAchievesXMinusEps = subExpr;
-                            } else {
-                                pointAchievesXMinusEps = pointAchievesXMinusEps && subExpr;
-                            }
-                        }
-                        context.smtSolver->add(!pointAchievesXMinusEps);
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::addNewSimplexPoint(FacetAnalysisContext& context, PointId const& pointId, bool performCheck) {
+                auto const& coordinates = pointset.getPoint(pointId).get();
+                storm::expressions::Expression pointAchievesXMinusEps;
+                for (uint64_t i = 0; i < coordinates.size(); ++i) {
+                    storm::expressions::Expression subExpr = context.xMinusEps[i] <= context.expressionManager->rational(coordinates[i]);
+                    if (i == 0) {
+                        pointAchievesXMinusEps = subExpr;
+                    } else {
+                        pointAchievesXMinusEps = pointAchievesXMinusEps && subExpr;
                     }
-                    
-                    context.smtSolver->push();
-                    
+                }
+                context.smtSolver->add(!pointAchievesXMinusEps);
+                if (performCheck) {
                     auto smtCheckResult = context.smtSolver->check();
                     if (smtCheckResult == storm::solver::SmtSolver::CheckResult::Unsat) {
-                        // There is no point x that violates the 'precision criterion' for this facet
+                        // For all points x, there is a cached point that dominates or is equal to (x-eps).
+                        // (we have a constraint pointAchievesXminusEps that does not not hold (double negation)
                         return true;
                     } else {
                         STORM_LOG_THROW(smtCheckResult == storm::solver::SmtSolver::CheckResult::Sat, storm::exceptions::UnexpectedException, "The smt solver did not yield sat or unsat.");
+                        // there is a point x such that (x-eps) is not dominated by or equal to a cached point.
                         return false;
                     }
+                } else {
+                    return false;
                 }
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::findAndCheckCachedPoints(Environment const& env, FacetAnalysisContext& context) {
+                /*
+                Polytope inducedPoly = context.facet.getInducedSimplex(pointset);
+                pointset.collectPointsInPolytope(context.collectedPoints, inducedPoly);
+                uint64_t numNewPoints = context.collectedPoints.size();
+                STORM_LOG_ASSERT(numNewPoints >= context.facet.getNumberOfPoints(), "Did not find all points on the facet");
+                // return true iff for all points x there is a cached point that dominates or is equal to (x-eps).
+                for (auto const& pId : context.collectedPoints) {
+                    --numNewPoints;
+                    if (numNewPoints == 0) {
+                        return addNewSimplexPoint(context, pId, true);
+                    } else {
+                        addNewSimplexPoint(context, pId, false);
+                    }
+                }
+                 */
+                STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Reached code that should be unreachable...");
+                
                 return false;
             }
             
             template <class SparseModelType, typename GeometryValueType>
             bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::analyzePointsOnFacet(Environment const& env, FacetAnalysisContext& context) {
-            
-                // TODO
+                // Enumerate all points on the facet by creating a sub-MDP
+                
+                // TODO: Enumerate them using the scheduler evaluator, ie create a class similar to the weight vector checker
+                
                 return false;
             }
             
             template <class SparseModelType, typename GeometryValueType>
             bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::analyzePointsInSimplex(Environment const& env, FacetAnalysisContext& context) {
-                // TODO
+                auto const& pointIds = context.facet.getPoints();
+                std::vector<typename DetSchedsSimplexChecker<SparseModelType, GeometryValueType>::Point> pointsOnFacet;
+                pointsOnFacet.reserve(pointIds.size());
+                for (auto const& pId : pointIds) {
+                    pointsOnFacet.push_back(pointset.getPoint(pId).get());
+                }
+                
+//                simplexChecker->setSimplex(context.facet.getInducedSimplex(pointset), context.facet.getHalfspace().normalVector(), pointsOnFacet);
+                //for (auto const& pointInSimplex : context.collectedPoints) {
+                //    simplexChecker->addAchievablePoint(pointset.getPoint(pointInSimplex).get());
+                //}
+                
+
                 return false;
             }
             

src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h

@@ -6,6 +6,7 @@
 #include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h"
 #include "storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h"
 #include "storm/modelchecker/multiobjective/deterministicScheds/DetSchedsWeightVectorChecker.h"
+#include "storm/modelchecker/multiobjective/deterministicScheds/DetSchedsSimplexChecker.h"
 
 #include "storm/storage/geometry/Polytope.h"
 #include "storm/storage/geometry/Halfspace.h"
@@ -114,10 +115,8 @@ namespace storm {
                     
                     /*!
                      * Creates a polytope that captures all points that lie 'under' the facet
-                     * More precisely, the vertices of the polytope are the points on the facet
-                     * and point p with p_i = min {x_i | x lies on facet}
                      */
-                    Polytope const& getInducedSimplex(Pointset const& pointset);
+                    Polytope const& getInducedSimplex(Pointset const& pointset, std::vector<GeometryValueType> const& referenceCoordinates);
                     
                     
 
@@ -160,7 +159,7 @@ namespace storm {
                 void addHalfspaceToOverApproximation(Environment const& env, std::vector<GeometryValueType> const& normalVector, Point const& pointOnHalfspace);
                 
                 /*!
-                 * Adds a polytope which consists of unachievable point
+                 * Adds a polytope which consists of unachievable points
                  */
                 void addUnachievableArea(Environment const& env, Polytope const& area);
                 
@@ -170,6 +169,11 @@ namespace storm {
                  */
                 void initializeFacets(Environment const& env);
                 
+                /*!
+                 *  Gets reference coordinates used to subdividing the downwardclosure
+                 */
+                std::vector<GeometryValueType> getReferenceCoordinates() const;
+                
                 /*!
                  * Checks the precision of the given Facet and returns true, if no further processing of the facet is necessary
                  */
@@ -199,6 +203,15 @@ namespace storm {
                  */
                 bool optimizeAndSplitFacet(Environment const& env, Facet& f);
                 
+                /*!
+                 * Adds a new point that lies within the induced simplex of the given facet to the analysis context.
+                 * @param context the analysis context
+                 * @param pointId the id of the given point.
+                 * @param performCheck if true, it is checked whether the facet is sufficiently precise now. If false, no check is performed.
+                 * @return true iff performCheck is true and the facet is sufficiently precise.
+                 */
+                bool addNewSimplexPoint(FacetAnalysisContext& context, PointId const& pointId, bool performCheck);
+                
                 /*!
                  * Finds all points that lie within the induced Simplex of the given facet.
                  * Returns true if the facet is sufficiently precise when considering all added points
@@ -209,9 +222,6 @@ namespace storm {
                  * Finds points that lie on the facet
                  * Returns true if the facet has been analyzed sufficiently precise.
                  * If false is returned, it means that *all* points that lie on the facet have been analyzed but the analysis is still not precise enough
-                 *
-                 * use smt to find an eps-box in the simplex that does not contain a point. Add new points until one in the box is found. repeat.
-                 * stop when no more points or boxes can be found.
                  */
                 bool analyzePointsOnFacet(Environment const& env, FacetAnalysisContext& context);
                 
@@ -224,6 +234,7 @@ namespace storm {
                 
                 std::shared_ptr<MultiObjectiveSchedulerEvaluator<SparseModelType>> schedulerEvaluator;
                 std::shared_ptr<DetSchedsWeightVectorChecker<SparseModelType>> weightVectorChecker;
+                std::shared_ptr<DetSchedsSimplexChecker<SparseModelType, GeometryValueType>> simplexChecker;
                 std::shared_ptr<SparseModelType> const& model;
                 uint64_t originalModelInitialState;
                 std::vector<Objective<ModelValueType>> const& objectives;

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

@@ -100,8 +100,8 @@ namespace storm {
                     }
                     
                     auto detModel = mdp ? mdp->applyScheduler(scheduler, false) : model.applyScheduler(scheduler, false);
-                    STORM_LOG_ASSERT(detModel->isOfType(storm::models::ModelType::Dtmc), "Model is of unexpected type.");
-                    auto const& dtmc = *(detModel->template as<storm::models::sparse::Dtmc<ValueType>>());
+                    detModel->getTransitionMatrix().makeRowGroupingTrivial();
+                    storm::models::sparse::Dtmc<ValueType> dtmc(std::move(detModel->getTransitionMatrix()), std::move(detModel->getStateLabeling()), std::move(detModel->getRewardModels()));
                     for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                         storm::modelchecker::CheckTask<storm::logic::Formula, ValueType> task(*this->objectives[objIndex].formula, false);
                         auto res = storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<ValueType>>(dtmc).check(env, task);
@@ -151,6 +151,13 @@ namespace storm {
                 return results[objIndex];
             }
             
+            
+            template <class ModelType>
+            typename MultiObjectiveSchedulerEvaluator<ModelType>::ValueType const& MultiObjectiveSchedulerEvaluator<ModelType>::getSchedulerIndependentStateResult(uint64_t objIndex, uint64_t state) const {
+                STORM_LOG_ASSERT(getSchedulerIndependentStates(objIndex).get(state), "Tried to get scheduler-independent result for a scheduler-dependent state.");
+                return results[objIndex][state];
+            }
+            
             template <class ModelType>
             std::vector<std::vector<typename MultiObjectiveSchedulerEvaluator<ModelType>::ValueType>> const& MultiObjectiveSchedulerEvaluator<ModelType>::getResults() const {
                 STORM_LOG_ASSERT(currSchedHasBeenChecked, "Tried to get results for a scheduler that has not yet been analyzed.");

src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h

@@ -31,6 +31,7 @@ namespace storm {
                 // Retrieve the results after calling check.
                 std::vector<std::vector<ValueType>> const& getResults() const;
                 std::vector<ValueType> const& getResultForObjective(uint64_t objIndex) const;
+                ValueType const& getSchedulerIndependentStateResult(uint64_t objIndex, uint64_t state) const;
                 storm::storage::BitVector const& getSchedulerIndependentStates(uint64_t objIndex) const;
                 std::vector<ValueType> getInitialStateResults() const;
                 ModelType const& getModel() const;