Storm
A Modern Probabilistic Model Checker
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SparsePcaaQuery.cpp
Go to the documentation of this file.
2
5#include "storm/io/export.h"
16
18
19namespace storm {
20namespace modelchecker {
21namespace multiobjective {
22
23template<class SparseModelType, typename GeometryValueType>
26 STORM_PRINT_AND_LOG("Pareto Curve Approximation Algorithm terminated after " << this->refinementSteps.size() << " refinement steps.\n");
27 }
28}
29
30template<class SparseModelType, typename GeometryValueType>
39
40template<class SparseModelType, typename GeometryValueType>
42 Point const& pointToBeSeparated) {
43 STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by "
44 << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(pointToBeSeparated)) << ".");
45
46 if (underApproximation->isEmpty()) {
47 // In this case, every weight vector is separating
48 uint_fast64_t objIndex = diracWeightVectorsToBeChecked.getNextSetIndex(0) % pointToBeSeparated.size();
49 WeightVector result(pointToBeSeparated.size(), storm::utility::zero<GeometryValueType>());
50 result[objIndex] = storm::utility::one<GeometryValueType>();
51 diracWeightVectorsToBeChecked.set(objIndex, false);
52 return result;
53 }
54
55 // Reaching this point means that the underApproximation contains halfspaces. The seperating vector has to be the normal vector of one of these halfspaces.
56 // We pick one with maximal distance to the given point. However, Dirac weight vectors that only assign a non-zero weight to a single objective take
57 // precedence.
58 STORM_LOG_ASSERT(!underApproximation->contains(pointToBeSeparated),
59 "Tried to find a separating point but the point is already contained in the underApproximation");
60 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = underApproximation->getHalfspaces();
61 uint_fast64_t farestHalfspaceIndex = halfspaces.size();
62 GeometryValueType farestDistance = -storm::utility::one<GeometryValueType>();
63 bool foundSeparatingDiracVector = false;
64 for (uint_fast64_t halfspaceIndex = 0; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) {
65 GeometryValueType distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated);
66 if (!storm::utility::isZero(distance)) {
67 storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero<GeometryValueType>(halfspaces[halfspaceIndex].normalVector());
68 bool isSingleObjectiveVector =
69 nonZeroVectorEntries.getNumberOfSetBits() == 1 && diracWeightVectorsToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0));
70 // Check if this halfspace is a better candidate than the current one
71 if ((!foundSeparatingDiracVector && isSingleObjectiveVector) ||
72 (foundSeparatingDiracVector == isSingleObjectiveVector && distance > farestDistance)) {
73 foundSeparatingDiracVector = foundSeparatingDiracVector || isSingleObjectiveVector;
74 farestHalfspaceIndex = halfspaceIndex;
75 farestDistance = distance;
76 }
77 }
78 }
79 if (foundSeparatingDiracVector) {
80 diracWeightVectorsToBeChecked &= storm::utility::vector::filterZero<GeometryValueType>(halfspaces[farestHalfspaceIndex].normalVector());
81 }
82
83 STORM_LOG_THROW(farestHalfspaceIndex < halfspaces.size(), storm::exceptions::UnexpectedException, "There is no seperating vector.");
84 STORM_LOG_DEBUG("Found separating weight vector: "
85 << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(halfspaces[farestHalfspaceIndex].normalVector()))
86 << ".");
87 return halfspaces[farestHalfspaceIndex].normalVector();
88}
89
90template<class SparseModelType, typename GeometryValueType>
92 // Normalize the direction vector so that the entries sum up to one
94 direction, storm::utility::one<GeometryValueType>() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero<GeometryValueType>()));
95 weightVectorChecker->check(env, storm::utility::vector::convertNumericVector<typename SparseModelType::ValueType>(direction));
96 STORM_LOG_DEBUG("weighted objectives checker result (under approximation) is " << storm::utility::vector::toString(
97 storm::utility::vector::convertNumericVector<double>(weightVectorChecker->getUnderApproximationOfInitialStateResults())));
98 RefinementStep step;
99 step.weightVector = direction;
100 step.lowerBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getUnderApproximationOfInitialStateResults());
101 step.upperBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getOverApproximationOfInitialStateResults());
102 // For the minimizing objectives, we need to scale the corresponding entries with -1 as we want to consider the downward closure
103 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
104 if (storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType())) {
105 step.lowerBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
106 step.upperBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
107 }
108 }
109 refinementSteps.push_back(std::move(step));
110
111 updateOverApproximation();
112 updateUnderApproximation();
113}
114
115template<class SparseModelType, typename GeometryValueType>
118 refinementSteps.back().weightVector, storm::utility::vector::dotProduct(refinementSteps.back().weightVector, refinementSteps.back().upperBoundPoint));
119
120 // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation,
121 // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case.
122 GeometryValueType maximumOffset = h.offset();
123 for (auto const& step : refinementSteps) {
124 maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), step.lowerBoundPoint));
125 }
126 if (maximumOffset > h.offset()) {
127 // We correct the issue by shifting the halfspace such that it contains the underapproximation
128 h.offset() = maximumOffset;
129 STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by "
130 << storm::utility::convertNumber<double>(h.invert().euclideanDistance(refinementSteps.back().upperBoundPoint)) << ".");
131 }
132 overApproximation = overApproximation->intersection(h);
133 STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true));
134}
135
136template<class SparseModelType, typename GeometryValueType>
138 std::vector<Point> paretoPoints;
139 paretoPoints.reserve(refinementSteps.size());
140 for (auto const& step : refinementSteps) {
141 paretoPoints.push_back(step.lowerBoundPoint);
142 }
144 STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true));
145}
146
147template<class SparseModelType, typename GeometryValueType>
149 return env.modelchecker().multi().isMaxStepsSet() && this->refinementSteps.size() >= env.modelchecker().multi().getMaxSteps();
150}
151
152template<typename SparseModelType, typename GeometryValueType>
154 STORM_LOG_ERROR_COND(objectives.size() == 2, "Exporting plot requested but this is only implemented for the two-dimensional case.");
155
156 auto transformedUnderApprox = transformObjectivePolytopeToOriginal(this->objectives, underApproximation);
157 auto transformedOverApprox = transformObjectivePolytopeToOriginal(this->objectives, overApproximation);
158
159 // Get pareto points as well as a hyperrectangle that is used to guarantee that the resulting polytopes are bounded.
161 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()),
162 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()));
163 std::vector<std::vector<GeometryValueType>> paretoPoints;
164 paretoPoints.reserve(refinementSteps.size());
165 for (auto const& step : refinementSteps) {
166 paretoPoints.push_back(transformObjectiveValuesToOriginal(this->objectives, step.lowerBoundPoint));
167 boundaries.enlarge(paretoPoints.back());
168 }
169 auto underApproxVertices = transformedUnderApprox->getVertices();
170 for (auto const& v : underApproxVertices) {
171 boundaries.enlarge(v);
172 }
173 auto overApproxVertices = transformedOverApprox->getVertices();
174 for (auto const& v : overApproxVertices) {
175 boundaries.enlarge(v);
176 }
177
178 // Further enlarge the boundaries a little
179 storm::utility::vector::scaleVectorInPlace(boundaries.lowerBounds(), GeometryValueType(15) / GeometryValueType(10));
180 storm::utility::vector::scaleVectorInPlace(boundaries.upperBounds(), GeometryValueType(15) / GeometryValueType(10));
181
182 auto boundariesAsPolytope = boundaries.asPolytope();
183 std::vector<std::string> columnHeaders = {"x", "y"};
184
185 std::vector<std::vector<double>> pointsForPlotting;
187 underApproxVertices = transformedUnderApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
188 pointsForPlotting.reserve(underApproxVertices.size());
189 for (auto const& v : underApproxVertices) {
190 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
191 }
192 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathUnderApproximation().get(), pointsForPlotting, columnHeaders);
193 }
194
196 pointsForPlotting.clear();
197 overApproxVertices = transformedOverApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
198 pointsForPlotting.reserve(overApproxVertices.size());
199 for (auto const& v : overApproxVertices) {
200 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
201 }
202 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathOverApproximation().get(), pointsForPlotting, columnHeaders);
203 }
204
206 pointsForPlotting.clear();
207 pointsForPlotting.reserve(paretoPoints.size());
208 for (auto const& v : paretoPoints) {
209 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
210 }
211 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathParetoPoints().get(), pointsForPlotting, columnHeaders);
212 }
213}
214
215#ifdef STORM_HAVE_CARL
216template class SparsePcaaQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>;
217template class SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
218
219template class SparsePcaaQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
221#endif
222} // namespace multiobjective
223} // namespace modelchecker
224} // namespace storm
ModelCheckerEnvironment & modelchecker()
MultiObjectiveModelCheckerEnvironment & multi()
void performRefinementStep(Environment const &env, WeightVector &&direction)
SparsePcaaQuery(preprocessing::SparseMultiObjectivePreprocessorResult< SparseModelType > &preprocessorResult)
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > overApproximation
WeightVector findSeparatingVector(Point const &pointToBeSeparated)
std::unique_ptr< PcaaWeightVectorChecker< SparseModelType > > weightVectorChecker
bool maxStepsPerformed(Environment const &env) const
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > underApproximation
void exportPlotOfCurrentApproximation(Environment const &env) const
std::vector< Objective< typename SparseModelType::ValueType > > objectives
static std::unique_ptr< PcaaWeightVectorChecker< ModelType > > create(preprocessing::SparseMultiObjectivePreprocessorResult< ModelType > const &preprocessorResult)
A bit vector that is internally represented as a vector of 64-bit values.
Definition BitVector.h:18
uint_fast64_t getNextSetIndex(uint_fast64_t startingIndex) const
Retrieves the index of the bit that is the next bit set to true in the bit vector.
uint_fast64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
Halfspace< ValueType > invert() const
Definition Halfspace.h:67
std::vector< ValueType > const & normalVector() const
Definition Halfspace.h:112
ValueType const & offset() const
Definition Halfspace.h:120
void enlarge(std::vector< ValueType > const &point)
std::vector< ValueType > const & lowerBounds() const
std::shared_ptr< Polytope< ValueType > > asPolytope() const
std::vector< ValueType > const & upperBounds() const
static std::shared_ptr< Polytope< ValueType > > createDownwardClosure(std::vector< Point > const &points)
Creates the downward closure of the given points (i.e., the set { x | ex.
Definition Polytope.cpp:43
static std::shared_ptr< Polytope< ValueType > > createUniversalPolytope()
Creates the universal polytope (i.e., the set R^n)
Definition Polytope.cpp:27
static std::shared_ptr< Polytope< ValueType > > createEmptyPolytope()
Creates the empty polytope (i.e., emptyset)
Definition Polytope.cpp:32
#define STORM_LOG_WARN(message)
Definition logging.h:30
#define STORM_LOG_DEBUG(message)
Definition logging.h:23
#define STORM_LOG_ASSERT(cond, message)
Definition macros.h:11
#define STORM_LOG_ERROR_COND(cond, message)
Definition macros.h:52
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
#define STORM_PRINT_AND_LOG(message)
Definition macros.h:68
std::vector< GeometryValueType > transformObjectiveValuesToOriginal(std::vector< Objective< ValueType > > objectives, std::vector< GeometryValueType > const &point)
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > transformObjectivePolytopeToOriginal(std::vector< Objective< ValueType > > objectives, std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > const &polytope)
SettingsType const & getModule()
Get module.
bool constexpr minimize(OptimizationDirection d)
T dotProduct(std::vector< T > const &firstOperand, std::vector< T > const &secondOperand)
Computes the dot product (aka scalar product) and returns the result.
Definition vector.h:517
std::string toString(std::vector< ValueType > const &vector)
Output vector as string.
Definition vector.h:1298
void scaleVectorInPlace(std::vector< ValueType1 > &target, ValueType2 const &factor)
Multiplies each element of the given vector with the given factor and writes the result into the vect...
Definition vector.h:491
bool isZero(ValueType const &a)
Definition constants.cpp:41
LabParser.cpp.
Definition cli.cpp:18