diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp index 5c15256e8..c7f69cf8f 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp @@ -1,11 +1,13 @@ #include -#include +#include + #include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h" +#include "storm/storage/geometry/coordinates.h" #include "storm/models/sparse/MarkovAutomaton.h" #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/StandardRewardModel.h" - +#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/utility/export.h" @@ -126,7 +128,7 @@ namespace storm { } template - boost::optional::PointId> DeterministicParetoExplorer::Pointset::addPoint(Point const& point) { + boost::optional::PointId> DeterministicParetoExplorer::Pointset::addPoint(Point&& point) { // Find dominated and dominating points auto pointsIt = points.begin(); @@ -138,6 +140,10 @@ namespace storm { ++pointsIt; case Point::DominanceResult::Dominates: // Found a point in the set that is dominated by the new point, so we erase it + if (pointsIt->second.isParetoOptimal()) { + STORM_LOG_WARN("Potential precision issues: Found a point that dominates another point which was flagged as pareto optimal. Distance of points is " << storm::storage::geometry::euclideanDistance(pointsIt->second.get(), point.get())); + point.setParetoOptimal(true); + } pointsIt = points.erase(pointsIt); case Point::DominanceResult::Dominated: // The new point is dominated by another point. @@ -150,7 +156,7 @@ namespace storm { } } - points.emplace_hint(points.end(), currId, point); + points.emplace_hint(points.end(), currId, std::move(point)); return currId++; } @@ -159,11 +165,31 @@ namespace storm { return points.at(id); } + template + typename DeterministicParetoExplorer::Pointset::iterator_type DeterministicParetoExplorer::Pointset::begin() const { + return points.begin(); + } + + template + typename DeterministicParetoExplorer::Pointset::iterator_type DeterministicParetoExplorer::Pointset::end() const { + return points.end(); + } + template uint64_t DeterministicParetoExplorer::Pointset::size() const { return points.size(); } + template + typename DeterministicParetoExplorer::Polytope DeterministicParetoExplorer::Pointset::downwardClosure() const { + std::vector> pointsAsVector; + pointsAsVector.reserve(size()); + for (auto const& p : points) { + pointsAsVector.push_back(p.second.get()); + } + return storm::storage::geometry::Polytope::createDownwardClosure(std::move(pointsAsVector)); + } + template void DeterministicParetoExplorer::Pointset::collectPointsInPolytope(std::set& collectedPoints, Polytope const& polytope) { for (auto const& p : points) { @@ -194,8 +220,14 @@ namespace storm { // Intentionally left empty } + template + storm::storage::geometry::Halfspace const& DeterministicParetoExplorer::Facet::getHalfspace() const { + return halfspace; + } + template void DeterministicParetoExplorer::Facet::addPoint(PointId const& pointId) { + inducedSimplex = nullptr; paretoPointsOnFacet.push_back(pointId); } @@ -211,37 +243,41 @@ namespace storm { template - typename DeterministicParetoExplorer::Polytope DeterministicParetoExplorer::Facet::getInducedSimplex(Pointset const& pointset) const { - std::vector> 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 minPoint = vertices.back(); - - for (++pointIdIt; pointIdIt != pointIdItE; ++pointIdIt) { + typename DeterministicParetoExplorer::Polytope const& DeterministicParetoExplorer::Facet::getInducedSimplex(Pointset const& pointset) { + if (!inducedSimplex) { + std::vector> 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()); - auto pIt = vertices.back().begin(); - for (auto& mi : minPoint) { - mi = std::min(mi, *pIt); - ++pIt; + std::vector 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; + } + } + 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"); + if (dimensionsForDownwardClosure.empty()) { + inducedSimplex = storm::storage::geometry::Polytope::create(vertices); + } else { + inducedSimplex = storm::storage::geometry::Polytope::createSelectiveDownwardClosure(vertices, dimensionsForDownwardClosure); } } - 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"); - if (dimensionsForDownwardClosure.empty()) { - return storm::storage::geometry::Polytope::create(vertices); - } else { - return storm::storage::geometry::Polytope::createSelectiveDownwardClosure(vertices, dimensionsForDownwardClosure); - } + return inducedSimplex; } template - DeterministicParetoExplorer::DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) { - // Intentionally left empty + DeterministicParetoExplorer::DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) : model(preprocessorResult.preprocessedModel), objectives(preprocessorResult.objectives) { + schedulerEvaluator = std::make_shared>(preprocessorResult); + weightVectorChecker = std::make_shared>(schedulerEvaluator); } template @@ -269,44 +305,183 @@ namespace storm { template void DeterministicParetoExplorer::initializeFacets(Environment const& env) { - // TODO + for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { + std::vector weightVector(objectives.size(), storm::utility::zero()); + if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { + weightVector[objIndex] = -storm::utility::one(); + } else { + weightVector[objIndex] = storm::utility::one(); + } + 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) { + if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { + (*pIt)[objIndex] *= -storm::utility::one(); + } + } + Point p(storm::utility::vector::convertNumericVector(*pIt)); + if (last) { + last = false; + p.setParetoOptimal(true); + // Adapt the overapproximation + std::vector normalVector(objectives.size(), storm::utility::zero()); + normalVector[objIndex] = storm::utility::one(); + GeometryValueType offset = storm::utility::vector::dotProduct(normalVector, p.get()); + storm::storage::geometry::Halfspace overApproxHalfspace(std::move(normalVector), std::move(offset)); + overApproximation = overApproximation->intersection(overApproxHalfspace); + } + pointset.addPoint(std::move(p)); + } + } + + auto initialHalfspaces = pointset.downwardClosure()->getHalfspaces(); + for (auto& h : initialHalfspaces) { + Facet f(std::move(h)); + for (auto const& p : pointset) { + if (f.getHalfspace().isPointOnBoundary(p.second.get())) { + f.addPoint(p.first); + } + } + STORM_LOG_ASSERT(std::count(f.getHalfspace().normalVector().begin(), f.getHalfspace().normalVector().end(), storm::utility::zero()) + f.getNumberOfPoints() == objectives.size(), "Unexpected number of points on facet."); + if (!checkFacetPrecision(env, f)) { + unprocessedFacets.push(std::move(f)); + } + } } template - bool DeterministicParetoExplorer::checkFacetPrecision(Environment const& env, Facet const& f) { - // TODO - return false; + bool DeterministicParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f) { + auto const& inducedSimplex = f.getInducedSimplex(pointset); + + GeometryValueType eps = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); + // get a polytope that contains exactly the points y, such that y+eps is in the induced simplex + std::vector offsetVector(objectives.size(), -eps); + auto shiftedSimplex = inducedSimplex->shift(offsetVector); + + // 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 - void DeterministicParetoExplorer::processFacet(Environment const& env, Facet const& f) { + bool DeterministicParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f, std::set const& collectedSimplexPoints) { // TODO - + return false; } template - bool DeterministicParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet const& f) { - // TODO - return false; + void DeterministicParetoExplorer::processFacet(Environment const& env, Facet& f) { + if (optimizeAndSplitFacet(env,f)) { + return; + } + std::set collectedPoints; + + if (findAndCheckCachedPoints(env, f, collectedPoints)) { + return; + } + + if (analyzePointsOnFacet(env, f, collectedPoints)) { + return; + } + + if (analyzePointsInSimplex(env, f, collectedPoints)) { + 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."); + } + template + bool DeterministicParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet& f) { + // Obtain the correct weight vector + auto weightVector = storm::utility::vector::convertNumericVector(f.getHalfspace().normalVector()); + for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { + weightVector[objIndex] *= -storm::utility::one(); + } + } + + // Invoke optimization and insert the explored points + boost::optional optPointId; + 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) { + if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { + (*pIt)[objIndex] *= -storm::utility::one(); + } + } + Point p(storm::utility::vector::convertNumericVector(*pIt)); + if (last) { + last = false; + p.setParetoOptimal(true); + // Adapt the overapproximation + GeometryValueType offset = storm::utility::vector::dotProduct(f.getHalfspace().normalVector(), p.get()); + storm::storage::geometry::Halfspace overApproxHalfspace(f.getHalfspace().normalVector(), std::move(offset)); + overApproximation = overApproximation->intersection(overApproxHalfspace); + optPointId = pointset.addPoint(std::move(p)); + } else { + pointset.addPoint(std::move(p)); + } + } + + // Potentially generate new facets + if (optPointId) { + auto const& optPoint = pointset.getPoint(*optPointId); + 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 + std::vector> vertices; + vertices.push_back(optPoint.get()); + for (auto const& pId : f.getPoints()) { + vertices.push_back(pointset.getPoint(pId).get()); + } + auto newHalfspaceCandidates = storm::storage::geometry::Polytope::createSelectiveDownwardClosure(vertices, storm::utility::vector::filterZero(f.getHalfspace().normalVector()))->getHalfspaces(); + for (auto& h : newHalfspaceCandidates) { + if (!storm::utility::vector::hasNegativeEntry(h.normalVector())) { + STORM_LOG_ASSERT(h.isPointOnBoundary(optPoint.get()), "Unexpected facet found while splitting."); + Facet fNew(std::move(h)); + fNew.addPoint(optPointId.get()); + auto vertexIt = vertices.begin(); + ++vertexIt; + for (auto const& pId : f.getPoints()) { + assert(pointset.getPoint(pId).get() == *vertexIt); + if (fNew.getHalfspace().isPointOnBoundary(*vertexIt)) { + fNew.addPoint(pId); + } + ++vertexIt; + } + assert(vertexIt == vertices.end()); + unprocessedFacets.push(std::move(fNew)); + } + } + return true; + } + } else { + // If the 'optimal point' was dominated by an existing point, we can not split the facet any further. + return false; + } } template - bool DeterministicParetoExplorer::findAndCheckCachedPoints(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints) { + bool DeterministicParetoExplorer::findAndCheckCachedPoints(Environment const& env, Facet& f, std::set& collectedPoints) { // TODO return false; } template - bool DeterministicParetoExplorer::analyzePointsOnFacet(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints) { + bool DeterministicParetoExplorer::analyzePointsOnFacet(Environment const& env, Facet& f, std::set& collectedPoints) { // TODO return false; } template - bool DeterministicParetoExplorer::analyzePointsInSimplex(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints) { - + bool DeterministicParetoExplorer::analyzePointsInSimplex(Environment const& env, Facet& f, std::set& collectedPoints) { // TODO return false; } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h index 3ca0b85d0..e0215b6ce 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h @@ -3,6 +3,9 @@ #include #include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DetSchedsWeightVectorChecker.h" + #include "storm/storage/geometry/Polytope.h" #include "storm/storage/geometry/Halfspace.h" #include "storm/modelchecker/results/CheckResult.h" @@ -19,6 +22,7 @@ namespace storm { public: typedef uint64_t PointId; typedef typename std::shared_ptr> Polytope; + typedef typename SparseModelType::ValueType ModelValueType; class Point { public: @@ -51,6 +55,10 @@ namespace storm { class Pointset { public: + + typedef typename std::map::const_iterator iterator_type; + + Pointset(); /*! @@ -60,18 +68,26 @@ namespace storm { * Erases all points in the set, that are dominated by the given point. * If the same point is already contained in the set, its id is returned */ - boost::optional addPoint(Point const& point); + boost::optional addPoint(Point&& point); /*! * Returns the point with the given ID */ Point const& getPoint(PointId const& id) const; + iterator_type begin() const; + iterator_type end() const; + /*! * Returns the number of points currently contained in the set */ uint64_t size() const; + /*! + * Returns the downward closure of the contained points. + */ + Polytope downwardClosure() const; + void collectPointsInPolytope(std::set& collectedPoints, Polytope const& polytope); void printToStream(std::ostream& out, bool includeIDs = true, bool convertToDouble = false); @@ -85,6 +101,7 @@ namespace storm { public: Facet(storm::storage::geometry::Halfspace const& halfspace); Facet(storm::storage::geometry::Halfspace&& halfspace); + storm::storage::geometry::Halfspace const& getHalfspace() const; void addPoint(PointId const& pointId); std::vector const& getPoints() const; uint64_t getNumberOfPoints() const; @@ -94,13 +111,14 @@ namespace storm { * 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 getInducedSimplex(Pointset const& pointset) const; + Polytope const& getInducedSimplex(Pointset const& pointset); private: storm::storage::geometry::Halfspace halfspace; std::vector paretoPointsOnFacet; + Polytope inducedSimplex; }; @@ -126,7 +144,13 @@ namespace storm { /*! * Checks the precision of the given Facet and returns true, if no further processing of the facet is necessary */ - bool checkFacetPrecision(Environment const& env, Facet const& f); + bool checkFacetPrecision(Environment const& env, Facet& f); + + /*! + * Checks the precision of the given Facet and returns true, if no further processing of the facet is necessary. + * Also takes the given points within the simplex of the facet into account + */ + bool checkFacetPrecision(Environment const& env, Facet& f, std::set const& collectedSimplexPoints); /*! Processes the given facet as follows: * 1. Optimize in the facet direction. Potentially, this adds new, unprocessed facets @@ -134,7 +158,7 @@ namespace storm { * 3. Find more points that lie on the facet * 4. Find all points that lie in the induced simplex or prove that there are none */ - void processFacet(Environment const& env, Facet const& f); + void processFacet(Environment const& env, Facet& f); /*! * Optimizes in the facet direction. If this results in a point that does not lie on the facet, @@ -142,13 +166,13 @@ 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 const& f); + bool optimizeAndSplitFacet(Environment const& env, Facet& f); /*! * 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 */ - bool findAndCheckCachedPoints(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints); + bool findAndCheckCachedPoints(Environment const& env, Facet& f, std::set& collectedPoints); /*! * Finds points that lie on the facet @@ -158,15 +182,19 @@ namespace storm { * 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, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints); + bool analyzePointsOnFacet(Environment const& env, Facet& f, std::set& collectedPoints); - bool analyzePointsInSimplex(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set& collectedPoints); + bool analyzePointsInSimplex(Environment const& env, Facet& f, std::set& collectedPoints); Pointset pointset; std::queue unprocessedFacets; Polytope overApproximation; std::vector unachievableAreas; + std::shared_ptr> schedulerEvaluator; + std::shared_ptr> weightVectorChecker; + std::shared_ptr const& model; + std::vector> const& objectives; }; }