diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsSimplexChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsSimplexChecker.h new file mode 100644 index 000000000..c8e7033aa --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsSimplexChecker.h @@ -0,0 +1,387 @@ +#pragma once + +#include + +#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 + class PolytopeTree { + typedef typename std::shared_ptr> Polytope; + typedef typename std::vector 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> 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(pointPlusEps); + for (auto const& coordinate : point) { + pointPlusEps.push_back(coordinate + eps); + } + auto downwardOfPoint = storm::storage::geometry::Polytope::createDownwardClosure({pointPlusEps}); + setMinus(downwardOfPoint); + } + + bool isEmpty() const { + return polytope == nullptr; + } + + void clear() { + children.clear(); + polytope = nullptr; + } + + Polytope getPolytope() const { + return polytope; + } + + std::vector& 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(vi) << ","; + } + s << "]\t"; + } + s << std::endl; + return s.str(); + } + + private: + + + + Polytope polytope; + std::vector> children; + }; + + template + class DetSchedsSimplexChecker { + public: + + typedef typename ModelType::ValueType ValueType; + typedef typename std::shared_ptr> Polytope; + typedef typename std::vector Point; + + DetSchedsSimplexChecker(std::shared_ptr> 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> check(storm::Environment const& env, std::vector const& weightVector, PolytopeTree& 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(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> checkRecursive(std::vector const& weightVector, PolytopeTree& 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 foundPoints; + std::vector 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(lpModel->getContinuousValue(objVar))); + } + auto halfspace = storm::storage::geometry::Halfspace(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::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 const& objective, uint64_t choiceIndex) { + auto const& model = schedulerEvaluator->getModel(); + storm::modelchecker::SparsePropositionalModelChecker mc(model); + auto const& formula = *objective.formula; + if (formula.isProbabilityOperatorFormula() && formula.getSubformula().isUntilFormula()) { + return storm::utility::zero(); + } 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()); + 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()); + 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("model"); + lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + initialStateResults.clear(); + + auto one = lpModel->getConstant(storm::utility::one()); + + // Create choice variables and assert that at least one choice is taken at each state. + std::vector 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 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 const& objective = objectives[objIndex]; + storm::storage::BitVector const& schedulerIndependentStates = schedulerEvaluator->getSchedulerIndependentStates(objIndex); + // Create state variables + std::vector 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> schedulerEvaluator; + + std::unique_ptr> lpModel; + std::vector initialStateResults; + std::vector currentObjectiveVariables; + + + storm::utility::Stopwatch swInit; + storm::utility::Stopwatch swCheck; + storm::utility::Stopwatch swLpSolve; + storm::utility::Stopwatch swLpBuild; + storm::utility::Stopwatch swAux; + }; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsWeightVectorChecker.cpp index 85aa0d212..bef009c4c 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DetSchedsWeightVectorChecker.cpp +++ b/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(); + // TODO: Choice rewards? for (auto const& entry : transitionMatrix.getRow(choiceOffset + choice)) { objValue += entry.getValue() * stateResults[objIndex][entry.getColumn()]; } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp index 337768578..024c5767e 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp +++ b/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 - typename DeterministicParetoExplorer::Polytope const& DeterministicParetoExplorer::Facet::getInducedSimplex(Pointset const& pointset) { + typename DeterministicParetoExplorer::Polytope const& DeterministicParetoExplorer::Facet::getInducedSimplex(Pointset const& pointset, std::vector const& referenceCoordinates) { 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()); - 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; - } + std::vector> 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>(preprocessorResult); weightVectorChecker = std::make_shared>(schedulerEvaluator); + simplexChecker = std::make_shared>(schedulerEvaluator); } template @@ -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(); } @@ -450,8 +437,22 @@ namespace storm { } } + template + std::vector DeterministicParetoExplorer::getReferenceCoordinates() const { + std::vector 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(value)); + } + return result; + } + + template bool DeterministicParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f) { + // TODO: + return false; + /* auto const& inducedSimplex = f.getInducedSimplex(pointset); GeometryValueType eps = storm::utility::convertNumber(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 bool DeterministicParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f, std::set const& collectedSimplexPoints) { - // TODO + assert(false); return false; } @@ -476,6 +478,26 @@ namespace storm { return; } + GeometryValueType eps = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); + eps += eps; // The unknown area (box) can actually have size 2*eps + PolytopeTree 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 typename DeterministicParetoExplorer::FacetAnalysisContext DeterministicParetoExplorer::createAnalysisContext(Environment const& env, Facet& f) { - FacetAnalysisContext res(f); + FacetAnalysisContext res(f); + /* res.expressionManager = std::make_shared(); 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(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(); } @@ -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(); } @@ -606,50 +627,80 @@ namespace storm { } template - bool DeterministicParetoExplorer::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::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 + bool DeterministicParetoExplorer::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 bool DeterministicParetoExplorer::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 bool DeterministicParetoExplorer::analyzePointsInSimplex(Environment const& env, FacetAnalysisContext& context) { - // TODO + auto const& pointIds = context.facet.getPoints(); + std::vector::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; } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h index 4bdb95fe2..42e01e52e 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h +++ b/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 const& referenceCoordinates); @@ -160,7 +159,7 @@ namespace storm { void addHalfspaceToOverApproximation(Environment const& env, std::vector 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 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> schedulerEvaluator; std::shared_ptr> weightVectorChecker; + std::shared_ptr> simplexChecker; std::shared_ptr const& model; uint64_t originalModelInitialState; std::vector> const& objectives; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.cpp index 99a1b5d79..f33bc7d5b 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.cpp +++ b/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>()); + detModel->getTransitionMatrix().makeRowGroupingTrivial(); + storm::models::sparse::Dtmc 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 task(*this->objectives[objIndex].formula, false); auto res = storm::modelchecker::SparseDtmcPrctlModelChecker>(dtmc).check(env, task); @@ -151,6 +151,13 @@ namespace storm { return results[objIndex]; } + + template + typename MultiObjectiveSchedulerEvaluator::ValueType const& MultiObjectiveSchedulerEvaluator::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 std::vector::ValueType>> const& MultiObjectiveSchedulerEvaluator::getResults() const { STORM_LOG_ASSERT(currSchedHasBeenChecked, "Tried to get results for a scheduler that has not yet been analyzed."); diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h b/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h index c2e4970aa..e3704d420 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h @@ -31,6 +31,7 @@ namespace storm { // Retrieve the results after calling check. std::vector> const& getResults() const; std::vector 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 getInitialStateResults() const; ModelType const& getModel() const;