started implementing deterministic scheduler finding approach for multi-objective model checking

src/storm/modelchecker/multiobjective/MultiObjectivePostprocessing.cpp

@@ -0,0 +1,73 @@
+#include "storm/modelchecker/multiobjective/MultiObjectivePostprocessing.h"
+
+#include "storm/adapters/RationalNumberAdapter.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template<typename ValueType, typename GeometryValueType>
+            std::vector<GeometryValueType> transformObjectiveValuesToOriginal(std::vector<Objective<ValueType>> objectives, std::vector<GeometryValueType> const& point) {
+                std::vector<GeometryValueType> result;
+                result.reserve(point.size());
+                for (uint_fast64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
+                    auto const& obj = objectives[objIndex];
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
+                        if (obj.considersComplementaryEvent) {
+                            result.push_back(storm::utility::one<GeometryValueType>() - point[objIndex]);
+                        } else {
+                            result.push_back(point[objIndex]);
+                        }
+                    } else {
+                        if (obj.considersComplementaryEvent) {
+                            result.push_back(storm::utility::one<GeometryValueType>() + point[objIndex]);
+                        } else {
+                            result.push_back(-point[objIndex]);
+                        }
+                    }
+                }
+                return result;
+            }
+            
+            template<typename ValueType, typename GeometryValueType>
+            std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> transformObjectivePolytopeToOriginal(std::vector<Objective<ValueType>> objectives, std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) {
+                if(polytope->isEmpty()) {
+                    return storm::storage::geometry::Polytope<GeometryValueType>::createEmptyPolytope();
+                }
+                if(polytope->isUniversal()) {
+                    return storm::storage::geometry::Polytope<GeometryValueType>::createUniversalPolytope();
+                }
+                uint_fast64_t numObjectives = objectives.size();
+                std::vector<std::vector<GeometryValueType>> transformationMatrix(numObjectives, std::vector<GeometryValueType>(numObjectives, storm::utility::zero<GeometryValueType>()));
+                std::vector<GeometryValueType> transformationVector;
+                transformationVector.reserve(numObjectives);
+                for(uint_fast64_t objIndex = 0; objIndex < numObjectives; ++objIndex) {
+                    auto const& obj = objectives[objIndex];
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
+                        if (obj.considersComplementaryEvent) {
+                            transformationMatrix[objIndex][objIndex] = -storm::utility::one<GeometryValueType>();
+                            transformationVector.push_back(storm::utility::one<GeometryValueType>());
+                        } else {
+                            transformationMatrix[objIndex][objIndex] = storm::utility::one<GeometryValueType>();
+                            transformationVector.push_back(storm::utility::zero<GeometryValueType>());
+                        }
+                    } else {
+                        if (obj.considersComplementaryEvent) {
+                            transformationMatrix[objIndex][objIndex] = storm::utility::one<GeometryValueType>();
+                            transformationVector.push_back(storm::utility::one<GeometryValueType>());
+                        } else {
+                            transformationMatrix[objIndex][objIndex] = -storm::utility::one<GeometryValueType>();
+                            transformationVector.push_back(storm::utility::zero<GeometryValueType>());
+                        }
+                    }
+                }
+                return polytope->affineTransformation(transformationMatrix, transformationVector);
+            }
+            
+            template std::vector<storm::RationalNumber> transformObjectiveValuesToOriginal(std::vector<Objective<double>> objectives, std::vector<storm::RationalNumber> const& point);
+            template std::shared_ptr<storm::storage::geometry::Polytope<storm::RationalNumber>> transformObjectivePolytopeToOriginal(std::vector<Objective<double>> objectives, std::shared_ptr<storm::storage::geometry::Polytope<storm::RationalNumber>> const& polytope);
+            template std::vector<storm::RationalNumber> transformObjectiveValuesToOriginal(std::vector<Objective<storm::RationalNumber>> objectives, std::vector<storm::RationalNumber> const& point);
+            template std::shared_ptr<storm::storage::geometry::Polytope<storm::RationalNumber>> transformObjectivePolytopeToOriginal(std::vector<Objective<storm::RationalNumber>> objectives, std::shared_ptr<storm::storage::geometry::Polytope<storm::RationalNumber>> const& polytope);
+        }
+    }
+}

src/storm/modelchecker/multiobjective/MultiObjectivePostprocessing.h

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <vector>
+#include <memory>
+
+#include "storm/modelchecker/multiobjective/Objective.h"
+#include "storm/storage/geometry/Polytope.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template<typename ValueType, typename GeometryValueType>
+            std::vector<GeometryValueType> transformObjectiveValuesToOriginal(std::vector<Objective<ValueType>> objectives, std::vector<GeometryValueType> const& point);
+            
+            template<typename ValueType, typename GeometryValueType>
+            std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> transformObjectivePolytopeToOriginal(std::vector<Objective<ValueType>> objectives, std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope);
+        
+        }
+    }
+}

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

@@ -0,0 +1,384 @@
+#include <sstream>
+#include <storm/models/sparse/MarkovAutomaton.h>
+
+#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h"
+#include "storm/models/sparse/MarkovAutomaton.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+
+
+#include "storm/utility/export.h"
+
+#include "storm/exceptions/InvalidOperationException.h"
+
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::Point(std::vector<GeometryValueType> const& coordinates) : coordinates(coordinates), paretoOptimal(false) {
+                STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported");
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::Point(std::vector<GeometryValueType>&& coordinates) : coordinates(std::move(coordinates)), paretoOptimal(false) {
+                STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported");
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            std::vector<GeometryValueType>& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::get() {
+                return coordinates;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            std::vector<GeometryValueType> const& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::get() const {
+                return coordinates;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            uint64_t DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::dimension() const {
+                STORM_LOG_ASSERT(!coordinates.empty(), "Points with dimension 0 are not supported");
+                return coordinates.size();
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::DominanceResult DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::getDominance(Point const& other) const {
+                STORM_LOG_ASSERT(this->dimension() == other.dimension(), "Non-Equal dimensions of points: [" << this->toString() << "] vs. [" << other.toString() << "]");
+                auto thisIt = this->get().begin();
+                auto otherIt = other.get().begin();
+                auto thisItE = this->get().end();
+                
+                // Find the first entry where the points differ
+                while (*thisIt == *otherIt) {
+                    ++thisIt;
+                    ++otherIt;
+                    if (thisIt == thisItE) {
+                        return DominanceResult::Equal;
+                    }
+                }
+                    
+                if (*thisIt > *otherIt) {
+                    // *this might dominate other
+                    for (++thisIt, ++otherIt; thisIt != thisItE; ++thisIt, ++otherIt) {
+                        if (*thisIt < *otherIt) {
+                            return DominanceResult::Incomparable;
+                        }
+                    }
+                    return DominanceResult::Dominates;
+                } else {
+                    assert(*thisIt < *otherIt);
+                    // *this might be dominated by other
+                    for (++thisIt, ++otherIt; thisIt != thisItE; ++thisIt, ++otherIt) {
+                        if (*thisIt > *otherIt) {
+                            return DominanceResult::Incomparable;
+                        }
+                    }
+                    return DominanceResult::Dominated;
+                }
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::setParetoOptimal(bool value) {
+                paretoOptimal = value;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::isParetoOptimal() const {
+                return paretoOptimal;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            std::string DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point::toString(bool convertToDouble) const {
+                std::stringstream out;
+                bool first = true;
+                for (auto const& pi : this->get()) {
+                    if (first) {
+                        first = false;
+                    } else {
+                        out << ", ";
+                    }
+                    if (convertToDouble) {
+                        out << pi;
+                    } else {
+                        out << pi;
+                    }
+                }
+                return out.str();
+            }
+            
+      //      template <class SparseModelType, typename GeometryValueType>
+      //      bool operator<(typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point const& lhs, typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point const& rhs) {
+      //          STORM_LOG_ASSERT(lhs.dimension() == rhs.dimension(), "Non-Equal dimensions of points: " << lhs << " vs. " << rhs);
+      //          for (uint64_t i = 0; i < lhs.dimension(); ++i) {
+      //              if (lhs.get()[i] < rhs.get()[i]) {
+      //                  return true;
+      //              } else if (lhs.get()[i] != rhs.get()[i]) {
+      //                  return false;
+      //              }
+      //          }
+      //          return false;
+      //      }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::Pointset() : currId(1) {
+                // Intentionally left empty
+            }
+        
+            template <class SparseModelType, typename GeometryValueType>
+            boost::optional<typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::PointId> DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::addPoint(Point const& point) {
+            
+                // 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:
+                            // Nothing to be done for this point
+                            ++pointsIt;
+                        case Point::DominanceResult::Dominates:
+                            // Found a point in the set that is dominated by the new point, so we erase it
+                            pointsIt = points.erase(pointsIt);
+                        case Point::DominanceResult::Dominated:
+                            // The new point is dominated by another point.
+                            return boost::none;
+                        case Point::DominanceResult::Equal:
+                            if (point.isParetoOptimal()) {
+                                pointsIt->second.setParetoOptimal();
+                            }
+                            return pointsIt->first;
+                    }
+                }
+                
+                points.emplace_hint(points.end(), currId, point);
+                return currId++;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Point const& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::getPoint(PointId const& id) const {
+                return points.at(id);
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            uint64_t DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::size() const {
+                return points.size();
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::collectPointsInPolytope(std::set<PointId>& collectedPoints, Polytope const& polytope) {
+                for (auto const& p : points) {
+                    if (polytope->contains(p.second.get())) {
+                        collectedPoints.insert(p.first);
+                    }
+                }
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Pointset::printToStream(std::ostream& out, bool includeIDs, bool convertToDouble) {
+                for (auto const& p : this->points) {
+                    if (includeIDs) {
+                        out << p.first << ": [" << p.second.toString(convertToDouble) << "]" << std::endl;
+                    } else {
+                        out << p.second.toString(convertToDouble) << std::endl;
+                    }
+                }
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::Facet(storm::storage::geometry::Halfspace<GeometryValueType> const& halfspace) : halfspace(halfspace) {
+                // Intentionally left empty
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::Facet(storm::storage::geometry::Halfspace<GeometryValueType>&& halfspace) : halfspace(std::move(halfspace)) {
+                // Intentionally left empty
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::addPoint(PointId const& pointId) {
+                paretoPointsOnFacet.push_back(pointId);
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            std::vector<typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::PointId> const& DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::getPoints() const {
+                return paretoPointsOnFacet;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            uint64_t DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::getNumberOfPoints() const {
+                return paretoPointsOnFacet.size();
+            }
+            
+            
+            template <class SparseModelType, typename GeometryValueType>
+            typename DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Polytope DeterministicParetoExplorer<SparseModelType, GeometryValueType>::Facet::getInducedSimplex(Pointset const& pointset) const {
+                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;
+                    }
+                }
+                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<GeometryValueType>::create(vertices);
+                } else {
+                    return storm::storage::geometry::Polytope<GeometryValueType>::createSelectiveDownwardClosure(vertices, dimensionsForDownwardClosure);
+                }
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            DeterministicParetoExplorer<SparseModelType, GeometryValueType>::DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult<SparseModelType>& preprocessorResult) {
+                // Intentionally left empty
+            }
+
+            template <class SparseModelType, typename GeometryValueType>
+            std::unique_ptr<CheckResult> DeterministicParetoExplorer<SparseModelType, GeometryValueType>::check(Environment const& env) {
+                
+                clean();
+                initializeFacets(env);
+                while (!unprocessedFacets.empty()) {
+                    Facet f = std::move(unprocessedFacets.front());
+                    unprocessedFacets.pop();
+                    processFacet(env, f);
+                }
+    
+                // Todo: Prepare check result
+                return nullptr;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::clean() {
+                pointset = Pointset();
+                unprocessedFacets = std::queue<Facet>();
+                overApproximation = storm::storage::geometry::Polytope<GeometryValueType>::createUniversalPolytope();
+                unachievableAreas.clear();
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::initializeFacets(Environment const& env) {
+                // TODO
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::checkFacetPrecision(Environment const& env, Facet const& f) {
+                // TODO
+                return false;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::processFacet(Environment const& env, Facet const& f) {
+                // TODO
+            
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet const& f) {
+                // TODO
+                return false;
+            
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::findAndCheckCachedPoints(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set<PointId>& collectedPoints) {
+                // TODO
+                return false;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::analyzePointsOnFacet(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set<PointId>& collectedPoints) {
+            
+                // TODO
+                return false;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            bool DeterministicParetoExplorer<SparseModelType, GeometryValueType>::analyzePointsInSimplex(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set<PointId>& collectedPoints) {
+            
+                // TODO
+                return false;
+            }
+            
+            template <class SparseModelType, typename GeometryValueType>
+            void DeterministicParetoExplorer<SparseModelType, GeometryValueType>::exportPlotOfCurrentApproximation(Environment const& env) {
+                /*
+                STORM_LOG_ERROR_COND(objectives.size()==2, "Exporting plot requested but this is only implemented for the two-dimensional case.");
+                
+                auto transformedUnderApprox = transformPolytopeToOriginalModel(underApproximation);
+                auto transformedOverApprox = transformPolytopeToOriginalModel(overApproximation);
+                
+                // Get pareto points as well as a hyperrectangle that is used to guarantee that the resulting polytopes are bounded.
+                storm::storage::geometry::Hyperrectangle<GeometryValueType> boundaries(std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()), std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()));
+                std::vector<std::vector<GeometryValueType>> paretoPoints;
+                paretoPoints.reserve(refinementSteps.size());
+                for(auto const& step : refinementSteps) {
+                    paretoPoints.push_back(transformPointToOriginalModel(step.lowerBoundPoint));
+                    boundaries.enlarge(paretoPoints.back());
+                }
+                auto underApproxVertices = transformedUnderApprox->getVertices();
+                for(auto const& v : underApproxVertices) {
+                    boundaries.enlarge(v);
+                }
+                auto overApproxVertices = transformedOverApprox->getVertices();
+                for(auto const& v : overApproxVertices) {
+                    boundaries.enlarge(v);
+                }
+                
+                //Further enlarge the boundaries a little
+                storm::utility::vector::scaleVectorInPlace(boundaries.lowerBounds(), GeometryValueType(15) / GeometryValueType(10));
+                storm::utility::vector::scaleVectorInPlace(boundaries.upperBounds(), GeometryValueType(15) / GeometryValueType(10));
+                
+                auto boundariesAsPolytope = boundaries.asPolytope();
+                std::vector<std::string> columnHeaders = {"x", "y"};
+                
+                std::vector<std::vector<double>> pointsForPlotting;
+                if (env.modelchecker().multi().getPlotPathUnderApproximation()) {
+                    underApproxVertices = transformedUnderApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
+                    pointsForPlotting.reserve(underApproxVertices.size());
+                    for(auto const& v : underApproxVertices) {
+                        pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
+                    }
+                    storm::utility::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathUnderApproximation().get(), pointsForPlotting, columnHeaders);
+                }
+                
+                if (env.modelchecker().multi().getPlotPathOverApproximation()) {
+                    pointsForPlotting.clear();
+                    overApproxVertices = transformedOverApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
+                    pointsForPlotting.reserve(overApproxVertices.size());
+                    for(auto const& v : overApproxVertices) {
+                        pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
+                    }
+                    storm::utility::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathOverApproximation().get(), pointsForPlotting, columnHeaders);
+                }
+                
+                if (env.modelchecker().multi().getPlotPathParetoPoints()) {
+                    pointsForPlotting.clear();
+                    pointsForPlotting.reserve(paretoPoints.size());
+                    for(auto const& v : paretoPoints) {
+                        pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
+                    }
+                    storm::utility::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathParetoPoints().get(), pointsForPlotting, columnHeaders);
+                }
+            };
+                 */
+            }
+            
+            template class DeterministicParetoExplorer<storm::models::sparse::Mdp<double>, storm::RationalNumber>;
+            template class DeterministicParetoExplorer<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
+            template class DeterministicParetoExplorer<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
+            template class DeterministicParetoExplorer<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
+        }
+    }
+}

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

@@ -0,0 +1,174 @@
+#pragma once
+
+#include <memory>
+
+#include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h"
+#include "storm/storage/geometry/Polytope.h"
+#include "storm/storage/geometry/Halfspace.h"
+#include "storm/modelchecker/results/CheckResult.h"
+
+namespace storm {
+    
+    class Environment;
+    
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template <class SparseModelType, typename GeometryValueType>
+            class DeterministicParetoExplorer {
+            public:
+                typedef uint64_t PointId;
+                typedef typename std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> Polytope;
+                
+                class Point {
+                public:
+                    Point(std::vector<GeometryValueType> const& coordinates);
+                    Point(std::vector<GeometryValueType>&& coordinates);
+                    
+                    std::vector<GeometryValueType> const& get() const;
+                    std::vector<GeometryValueType>& get();
+                    
+                    uint64_t dimension() const;
+                    
+                    enum class DominanceResult {
+                        Incomparable,
+                        Dominates,
+                        Dominated,
+                        Equal
+                    };
+                    DominanceResult getDominance(Point const& other) const;
+                    
+                    void setParetoOptimal(bool value = true);
+                    bool isParetoOptimal() const;
+                    
+                    std::string toString(bool convertToDouble = false) const;
+
+                private:
+                    std::vector<GeometryValueType> coordinates;
+                    bool paretoOptimal;
+                };
+                
+                
+                class Pointset {
+                public:
+                    Pointset();
+                    
+                    /*!
+                     * If the given point is not dominated by another point in the set, it is added
+                     * to the set and its ID is returned.
+                     * If the point is dominated by another point, boost::none is returned.
+                     * 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<PointId> addPoint(Point const& point);
+                    
+                    /*!
+                     * Returns the point with the given ID
+                     */
+                    Point const& getPoint(PointId const& id) const;
+                    
+                    /*!
+                     * Returns the number of points currently contained in the set
+                     */
+                    uint64_t size() const;
+                    
+                    void collectPointsInPolytope(std::set<PointId>& collectedPoints, Polytope const& polytope);
+                    
+                    void printToStream(std::ostream& out, bool includeIDs = true, bool convertToDouble = false);
+                    
+                private:
+                    std::map<PointId, Point> points;
+                    PointId currId;
+                };
+                
+                class Facet {
+                public:
+                    Facet(storm::storage::geometry::Halfspace<GeometryValueType> const& halfspace);
+                    Facet(storm::storage::geometry::Halfspace<GeometryValueType>&& halfspace);
+                    void addPoint(PointId const& pointId);
+                    std::vector<PointId> const& getPoints() const;
+                    uint64_t getNumberOfPoints() const;
+                    
+                    /*!
+                     * 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 getInducedSimplex(Pointset const& pointset) const;
+                    
+                    
+
+                private:
+                    storm::storage::geometry::Halfspace<GeometryValueType> halfspace;
+                    std::vector<PointId> paretoPointsOnFacet;
+                };
+                
+                
+                DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult<SparseModelType>& preprocessorResult);
+
+                virtual std::unique_ptr<CheckResult> check(Environment const& env);
+
+                void exportPlotOfCurrentApproximation(Environment const& env);
+                
+            private:
+                
+                /*!
+                 * Cleans up all cached results from a previous check call
+                 */
+                void clean();
+                
+                /*!
+                 *   Builds the initial facets by optimizing the objectives individually.
+                 *   Adds the facets that need further processing to unprocessedFacets
+                 */
+                void initializeFacets(Environment const& env);
+                
+                /*!
+                 * 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);
+                
+                /*! Processes the given facet as follows:
+                 * 1. Optimize in the facet direction. Potentially, this adds new, unprocessed facets
+                 * 2. Find points that have already been collected so far such that  they lie in the induced simplex of the facet.
+                 * 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);
+                
+                /*!
+                 * Optimizes in the facet direction. If this results in a point that does not lie on the facet,
+                 * 1. The new Pareto optimal point is added
+                 * 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);
+                
+                /*!
+                 * 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<PointId>& collectedPoints);
+                
+                /*!
+                 * 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, Facet const& f, Polytope const& inducedSimplex, std::set<PointId>& collectedPoints);
+                
+                bool analyzePointsInSimplex(Environment const& env, Facet const& f, Polytope const& inducedSimplex, std::set<PointId>& collectedPoints);
+                
+                Pointset pointset;
+                std::queue<Facet> unprocessedFacets;
+                Polytope overApproximation;
+                std::vector<Polytope> unachievableAreas;
+                
+            };
+            
+        }
+    }
+}

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

@@ -0,0 +1,32 @@
+#include "storm/modelchecker/multiobjective/deterministicScheds/MultiObjectiveSchedulerEvaluator.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template <class ModelType>
+            MultiObjectiveSchedulerEvaluator<ModelType>::MultiObjectiveSchedulerEvaluator(preprocessing::SparseMultiObjectivePreprocessorResult<ModelType>& preprocessorResult) {
+                // TODO
+            }
+            
+            template <class ModelType>
+            void MultiObjectiveSchedulerEvaluator<ModelType>::check(std::vector<uint64_t> const& scheduler) {
+                // TODO
+            }
+ 
+            template <class ModelType>
+            std::vector<typename MultiObjectiveSchedulerEvaluator<ModelType>::ValueType> const& MultiObjectiveSchedulerEvaluator<ModelType>::getResultForObjective(uint64_t objIndex) const {
+                return results[objIndex];
+            }
+            
+            template <class ModelType>
+            std::vector<ValueType> MultiObjectiveSchedulerEvaluator<ModelType>::getInitialStateResults() const {
+                std::vector<ValueType> res;
+                for (auto objResult : results) {
+                    res.push_back(objResult[initialState]);
+                }
+                return res;
+            }
+        }
+    }
+}

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

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <vector>
+
+#include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace multiobjective {
+            
+            template <class ModelType>
+            class MultiObjectiveSchedulerEvaluator {
+            public:
+                
+                typedef typename ModelType::ValueType ValueType;
+                
+                MultiObjectiveSchedulerEvaluator(preprocessing::SparseMultiObjectivePreprocessorResult<ModelType>& preprocessorResult);
+
+                void check(std::vector<uint64_t> const& scheduler);
+                
+                std::vector<ValueType> const& getResultForObjective(uint64_t objIndex) const;
+                std::vector<ValueType> getInitialStateResults() const;
+                
+            private:
+                std::vector<std::vector<ValueType>> results;
+                uint64_t initialState;
+            };
+        }
+    }
+}