diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp index c7f69cf8f..f3d2555c9 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp @@ -7,6 +7,8 @@ #include "storm/models/sparse/MarkovAutomaton.h" #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/StandardRewardModel.h" +#include "storm/modelchecker/multiobjective/MultiObjectivePostprocessing.h" +#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h" #include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/utility/export.h" @@ -101,7 +103,7 @@ namespace storm { out << ", "; } if (convertToDouble) { - out << pi; + out << storm::utility::convertNumber(pi); } else { out << pi; } @@ -128,7 +130,7 @@ namespace storm { } template - boost::optional::PointId> DeterministicParetoExplorer::Pointset::addPoint(Point&& point) { + boost::optional::PointId> DeterministicParetoExplorer::Pointset::addPoint(Environment const& env, Point&& point) { // Find dominated and dominating points auto pointsIt = points.begin(); @@ -138,13 +140,15 @@ namespace storm { case Point::DominanceResult::Incomparable: // Nothing to be done for this point ++pointsIt; + break; 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())); + STORM_LOG_WARN("Potential precision issues: Found a point that dominates another point which was flagged as pareto optimal. Distance of points is " << std::sqrt(storm::utility::convertNumber(storm::storage::geometry::squaredEuclideanDistance(pointsIt->second.get(), point.get())))); point.setParetoOptimal(true); } pointsIt = points.erase(pointsIt); + break; case Point::DominanceResult::Dominated: // The new point is dominated by another point. return boost::none; @@ -156,6 +160,10 @@ namespace storm { } } + if (env.modelchecker().multi().isPrintResultsSet()) { + std::cout << "## achievable point: [" << point.toString(true) << "]" << std::endl; + } + points.emplace_hint(points.end(), currId, std::move(point)); return currId++; } @@ -276,6 +284,7 @@ namespace storm { template DeterministicParetoExplorer::DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) : model(preprocessorResult.preprocessedModel), objectives(preprocessorResult.objectives) { + originalModelInitialState = *preprocessorResult.originalModel.getInitialStates().begin(); schedulerEvaluator = std::make_shared>(preprocessorResult); weightVectorChecker = std::make_shared>(schedulerEvaluator); } @@ -290,9 +299,14 @@ namespace storm { unprocessedFacets.pop(); processFacet(env, f); } - - // Todo: Prepare check result - return nullptr; + + std::vector>paretoPoints; + paretoPoints.reserve(pointset.size()); + for (auto const& p : pointset) { + paretoPoints.push_back(storm::utility::vector::convertNumericVector(transformObjectiveValuesToOriginal(objectives, p.second.get()))); + } + return std::make_unique>(originalModelInitialState, std::move(paretoPoints), + nullptr, nullptr); } template @@ -303,6 +317,60 @@ namespace storm { unachievableAreas.clear(); } + template + void DeterministicParetoExplorer::addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace) { + GeometryValueType offset = storm::utility::vector::dotProduct(normalVector, pointOnHalfspace.get()); + if (env.modelchecker().multi().isPrintResultsSet()) { + std::cout << "## unachievable halfspace: ["; + bool first = true; + for (auto const& xi : normalVector) { + if (first) { + first = false; + } else { + std::cout << ","; + } + std::cout << storm::utility::convertNumber(xi); + } + std::cout << "];[" << storm::utility::convertNumber(offset) << "]" << std::endl; + } + storm::storage::geometry::Halfspace overApproxHalfspace(normalVector, offset); + overApproximation = overApproximation->intersection(overApproxHalfspace); + } + + template + void DeterministicParetoExplorer::addUnachievableArea(Environment const& env, Polytope const& area) { + if (env.modelchecker().multi().isPrintResultsSet()) { + std::vector> vertices; + if (objectives.size() == 2) { + vertices = area->getVerticesInClockwiseOrder(); + } else { + vertices = area->getVertices(); + } + std::cout << "## unachievable polytope: "; + bool firstVertex = true; + for (auto const& v : vertices) { + if (firstVertex) { + firstVertex = false; + } else { + std::cout << ";"; + } + std::cout << "["; + bool firstEntry = true; + for (auto const& vi : v) { + if (firstEntry) { + firstEntry = false; + } else { + std::cout << ","; + } + std::cout << storm::utility::convertNumber(vi); + } + std::cout << "]"; + } + std::cout << std::endl; + } + unachievableAreas.push_back(area); + } + template void DeterministicParetoExplorer::initializeFacets(Environment const& env) { for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { @@ -327,11 +395,9 @@ namespace storm { // 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); + addHalfspaceToOverApproximation(env, normalVector, p); } - pointset.addPoint(std::move(p)); + pointset.addPoint(env, std::move(p)); } } @@ -390,7 +456,7 @@ namespace storm { } // 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 @@ -417,13 +483,10 @@ namespace storm { 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)); + addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), p); + optPointId = pointset.addPoint(env, std::move(p)); } else { - pointset.addPoint(std::move(p)); + pointset.addPoint(env, std::move(p)); } } @@ -456,7 +519,9 @@ namespace storm { ++vertexIt; } assert(vertexIt == vertices.end()); - unprocessedFacets.push(std::move(fNew)); + if (!checkFacetPrecision(env, fNew)) { + unprocessedFacets.push(std::move(fNew)); + } } } return true; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h index e0215b6ce..aeafc244c 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h @@ -68,7 +68,7 @@ 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&& point); + boost::optional addPoint(Environment const& env, Point&& point); /*! * Returns the point with the given ID @@ -135,6 +135,16 @@ namespace storm { */ void clean(); + /*! + * Intersects the overapproximation with the given halfspace + */ + void addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace); + + /*! + * Adds a polytope which consists of unachievable point + */ + void addUnachievableArea(Environment const& env, Polytope const& area); + /*! * Builds the initial facets by optimizing the objectives individually. * Adds the facets that need further processing to unprocessedFacets @@ -194,6 +204,7 @@ namespace storm { std::shared_ptr> schedulerEvaluator; std::shared_ptr> weightVectorChecker; std::shared_ptr const& model; + uint64_t originalModelInitialState; std::vector> const& objectives; }; diff --git a/src/storm/storage/geometry/coordinates.h b/src/storm/storage/geometry/coordinates.h index f2795440f..954df5497 100644 --- a/src/storm/storage/geometry/coordinates.h +++ b/src/storm/storage/geometry/coordinates.h @@ -6,7 +6,7 @@ namespace storm { namespace geometry { template - T euclideanDistance(std::vector const& p, std::vector const& q) { + T squaredEuclideanDistance(std::vector const& p, std::vector const& q) { STORM_LOG_ASSERT(p.size() == q.size(), "Invalid dimensions of input vectors."); T squaredSum = storm::utility::zero(); auto pIt = p.begin(); @@ -16,7 +16,7 @@ namespace storm { T diff = *pIt - *qIt; squaredSum += diff * diff; } - return storm::utility::sqrt(squaredSum); + return squaredSum; } }